Бүлэг 2 - ийн тайлбар

Тэгэхээр гол тэгшитгэл бол 2.5 бөгөөд энэ нь нэгдүгээр эрэмбийн бүрэн дифференциал тэгшитгэл юм. Энэ тэгшитгэлийг төгсгөлөг ялгаварын (finite difference method) аргаар бодъё.

Дээрх 2.9 тэгшитгэл нь програмчлалд бичигдсэн болно.

 гэдэг бол rh юм.
ce – г бид тэгшитгэл 2.4 – ээс тодорхойлох боломжтой. Мөн Т – г тэгшитгэл 2.5 – ийг ашиглан тодорхойлно. Бохирдлын концентрацын t = 0 нөхцөл дахь анхны утга  өгөгдвөл дараачын  хугацаанд концентраци  ямар байхыг тодорхойлох боломжтой болж байна.
Програмын бичлэг дээр асуух зүйл байвал тайлбарлах болно.

Бүлэг 3

Хайрцаган загварчлал дахь тооцоолон бодох аргууд /Numerical Solution for Box Model/

Яагаад хайрцаган загвар гэж нэрлэсэнийг нь мэдэхгүй юм. Мэдвэл тайлбар бичих болно. Эсвэл орчуулалгүй Бокс загвар /цаашид ХЗ гэе/ гэж ярьж болох л юм.

3.1 Ухагдахуун

ХЗ – ийн дифференциал тэгшитгэлийн энгийн тохиолдлыг нуурын усны чанар гэсэн сэдвээр авч үзэж тайлбарласан. Хэрэв оролтын эсвэл бохирын зарцуулга ямар нэгэн байдлаар өөрчлөгдөж байх нь ая холбогдол багатай юм. Түүгээрч барахгүй, ХЗ-ын маш олон хэрэглээнд сайхан шугаман тэгшитгэл тэр болгон тохиолдохгүй. Ерөнхийдөө, тэгшитгэлийг бодох тооцоолон бодох аргууд хэрэгтэй ба нуурын усны чанар гэсэн хэсэгт цухас дурьдагдсан болно.

C_e болон Т – г мэдэгдэж байгаа гэж санавал хугацааны функц нь тодорхой болно. Гэхдээ л хялбарчилахын тулд зарим нэг хэмжигдэхүүнийг тогтмолоор авъя. Хэрэв t_j хугацаан дахь концентаци c_j мэдэгдэж байвал тэгшитгэл (2.4) нь с – ийн өөрчлөлтийн утгыг өгнө.

Маш бага delta t хугацааны агшинд энэ утга тогтмол хэвээр үлдэнэ гэж авч үзвэл, тухайн интервалын дараах концентрацыг тодорхойлж болно.

t=to хугацаанд анхны нөхцөл co концентраци өгөгдөж эхлэхэд концентрацын шинж чанар нь delta t хугацаанд тооцоологдоно. Хугацааны алхам, эсвэл хугацааны интервал бага байх тусмаа илүү нарийвчлалтай байна гэдэг нь ойлгомжтой болж байна. Зураг 3.1 – д үүнийг дүрслэн харуулав. Эндээс урган гарч байгаа нэг асуулт нь юу вэ гэхээр хугацааны алхам яг ямар хэмжээтэй байж жинхэнэ нарийн үр дүнг авч болох вэ гэсэн асуулт юм.

! Computational Hydraulics

! Chapter 3 (Box Model - section 1)

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

    dt=3600       ! [sec] time interval

    nt=500         ! time step

   Q=1.0           ! [m3/sec] Qin, Qout

    V=6.e+5       ! [m3]

    Tr=3.3e+5    ! [sec] = 3.3day

    ci=5.0          ! [mg/l]=[g/m3]

    ds=100.0      ! [g/sec] L discharge

    ce=(ds+Q*ci)/(Q+V/Tr)        ! eq(2.4)

    T=V/(Q+V/Tr)                      ! eq(2.6)

    write(*,*) ce,T/(24.*3600) !; pause

    c0=10.0       ![mg/l]

     a=(ce-c0)/T

    c1=c0                                   ! initial concentration

    do n=1,nt                             ! advance of time

    time=n*dt                            ! time

            rh=(ce-c1)*(dt/T)         ! right hand side of eq (3.1).

    c2=rh+c1                             ! temporal integration

    c1=c2                                   ! next time step

    c3=a*time+c0                      !

    write(*,'(3f9.3)') time/T,c1/ce,c3/ce! output

    write(10,'(3f9.3)') time/T,c1/ce,c3/ce! output

    if(time/T.gt.3.0) goto 999

      enddo

      close(10)

  999 stop

  End program

Figure 1. Загварчлалын үр дүн /номны зурагтай ижил байна/

Эндээс dt гэдэг тоог dt = T*0.5(эсвэл 0.3 мөн 0.03) гэх мэтчилэн өгөгдлийг өгч тооцсон болно. Дээрх програмчлалын бичиглэлд dt=3600 гэж өгсөн байгаа.

Дээрх асуултын хариултыг дараачийн хэсгээс ойлгож магадгүй юм.

2. Тогтворжилт ба нарийвчлал

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

Хэрэвээ та энэ тэгшитгэлийг өмнөх бүлгийн 2.8 тэгшитгэлтэй харьцуулж үзвэл шал өөр мэт харагдана. Тэр тэгшитгэл нь энэ байна.    

Энэ байдал нь энгийн бөгөөд хэлбэр дүрсийн хувьд анхаарч үзэх нээх шаардлагагүй. Харин бүүр сайн ойртуулж харвал хоёулаа тодорхой нөхцөл дор нэг нэгэнтэйгээ ойр байх боломжтой харагдах ба гарцаагүй та ямар нөхцөл юм болдоо гэж бодож байгаа байх. Харваас хэмжих нэгжгүй параметр delta t/T нь тооцоолох аргын шинж чанараар тодорхойлогдоно. Энэ параметрын дөрөвхөн л тохиолдол байна. Үүнд:

Figure 2. Шийдиййн гуйвалт ба тогтворгүй байдал

Энэ нь тогтвортой байх нөхцөл нь ямар байх вэ гэдэгт анализ хийх боломж олгож байна.

Гэвч энэ тооцоолол нарийвчлэлтэй гэдгийг илэрхийлж чадахгүй. Хамгийн эцэст нь delta t/T байж болохоор хамгийн бага байх хэрэгтэй, гэвч хэрхэн бага байх вэ гэдгийг тоон болон аналитик шийдлийг хооронд нь харьцуулж байж л олох боломжтой. Жишээлбэл, “тооцооллын тайвшрах хугацаа” (numerical relaxation time) нь j*delta t хугацаанд дараах байдлаар тодорхойлогдоно.

Нарийчлалтай шийдын хувьд энэ нь нэгж утгатай ойролцоо байх хэрэгтэй. Тэгшитгэл 3.4 нь хугацааны алхамыг сонгох, нарийвчлалын критерийг тодорхой байдлаар урьдчилан танилцах боломжыг олгодог.

Хэд хэдэн хүчин зүйл хугацааны алхамд нөлөөлдөг. Энэ дотоод хугацааны хэмжээст харьцуулахад ч тэр маш бага байх хэрэгтэй ба ялгаагүй гадааг хугацааны хэмжээс Те – ийн хувьд ч тэр үүний 5 – 10% нь л байх хэрэгтэй. Маш хурдан өөрчлөлттэй тохиолдолд, энэ нь 3.4 тэгшитгэлээс гарах үр дүнтэй харьцуулахад хүндэвтэр шаардлага болно. Харин удаан өөрчлөлттэй үед, тэгшитгэл 3.4 нь илүү хязгаарлагдах ба маш нарийвчлалыг бий болгох тайвшрах хугацаа маань чухал биш болж харьцангуй том алдаанууд бий болж болзошгүй. Ямарч тохиолдолд, тогтворжилтын шаардлага нь судлагдсан байх хэрэгтэй.  

Өмнөх анализ нь тооцоолон бодсон шийд дахь алдааг авч үзсэн ба энд орхигдуулсан зүйлсийг илүү хэцүү тэгшитгэлийг ижил замаар тооцоолох үед ойлгох болно. Харамсалтай нь энэ бол анализ гэх мэтийг гүйцэтгэхэд үргэлж боломжтой байдаггүй. Шугаман биш тэгшитгэлд ч үргэлж боломжтой байх хувилбарууд байдаг. Гол нь ямар тэгшитгэл шийдвэрлэж байгаагаа л сайн мэдэж түүнд тохирсон аргыг л хэрэлгэх хэрэгтэй.

Энэ шинжилгээнд шийд (үл мэдэгдэгч) нь t_j хугацаан дахь төгсгөлөг ялгавар хэлбэртэй тэгшитгэлийн “суурь цэг” – ийг хамаарсан Тэйлорын цуваанаас тодорхойлогдоно.

Энэ бол “хасалтын алдаа” (truncation error) гэж нэрлэгдэх сүүлчийн илэрхийллээс бусад нь оргиналь дифференциал тэгшитгэл болно. Энэ тохиолдол нь delta t дахь нэгдүгээр эрэмбийн тэгшитгэл ба хэрэв delta t – г 1/2 хэмжээгээр багасгавал ижил хэмжээгээр алдаа нь багасах болно. Тэгшитгэл (3.1) – ийг хэрэглэж хасалтын алдааг тодорхойлж болно.

Тэгшитгэлийг (ce-c)/T – д харьцуулбал тэгшитгэл (3.4) – г тооцоолсонтой ижил харьцангуй алдаа гэх 1/2(delta t/T) – г өгнө. Хасалтын алдаа нь шийдэл хэдий тогтвортой болсон ч гэсэн алга болохгүй гэдгийг тэмдэглье. Тэгээд хугацааны алхамын утга нээх их чухал биш болно. Энэ нь тогтворжсон нөхцөлд тайвшрах хугацааг бий болгохоор шаардлагатай урт биш учраас өмнөх шинжилгээтэй тохирч байна.

3.3 Жишээ

Нуурын тухай авч үзсэн 2.1 хэсэгт БХХ- ийн хамгийн их хэмжээ бол нэг өдрийн цутгалт байсан. Тохирч буй тэнцвэрт концентраци нь 100мг/л байна. Хамгийн их концентацийг 50мг/л байхад зөвшөөрнө гэвэл дээрх хэмжээ хэр илүү гарах бол? Тайвшралтын хугацаа нь 2.5 өдөр гэдгийг тооцооллоор гаргаж ирсэн. Хугацааны алхам нь үр дүн бодитой нарийвчлалтай байхаар сонгогдох хэрэгтэй.

Гадаад хугацааны хэмжээс 1 өдөр (гэнэтийн өөрчлөлтийн хариу үйлдэл тайвшралтын хугацаанаас хамаарна). Нэг өдрийн дараа – хамгийн критик хугацаа – концентраци нарийвчлалтайгаар тодорхойлогдох хэрэгтэй, ойролцоогоор 1% - ийн зөрөөтэй, өөрөөр хэлбэл (1-delta t/T)^j  нь t=j*delta t=1 үед exp(-t/T) – ээс 1% - иас илүүгүй байх хэрэгтэй. Хэрэв delta t/T – ээр жаахан оролдвол 0.05 – ийн утганд хангалттай байхаар олдоно. Тэгэхээр delta t=3цаг байна. Энэ нь Te – тэй хамааруулж үзвэл боломжийн бага байна.

Симуляцийг хийхэд анхны нөхцөл шаардлагатай байна. Иймд co=30мг/л гэж авъя. Тооцоологдсон концентрацийг зураг 3.3 – д харуулав.

Критик утга нь 50мг/л – ээс илүү гарсан харагдаж байна. Энэ нь анхны нөхцөлөөс бага зэрэг хамаарсан байна.

3.4 Далд арга

3.1 хэсэгт байгаа Эйлэрийн арга нь хэрэв шийдэл нилээн аажим өөрчлөгдвөл (3.3) дахь тогтворжилтын критерээр хязгаарлагдана. Энэ бол энгийн ил арга (explicit method) юм. Энэ нь0<thetha<1  байх tj+thetha*delta t хугацаан дахь дифференциал тэгшитгэлийн ойролцоолоор хүчингүй болно. Энэ хугацаанд

Үл мэдэгдэгч cj+1 нь тэгшитгэлийн хоёр талд агуулагдаж байгаа энэ хэлбэрийн аргыг далд (implicit method) арга гэж нэрлэдэг. Алгебрын тэгшитгэл нь хугацааны алхам бүрт шийдвэрлэгдэх ёстой. Энэ тохиолдолд ганц шугаман тэгшитгэл байх хэрэгтэй боловч бусад тохиолдолд бэрх байж магадгүй. 3.5 хэсэгт дурьдагдсан жишээнд хасалтын алдааг дараах байдлаар гаргаж авсан.

Энэ нь нэгдүгээр эрэмбийн тэгшитгэл байсаар байгаа боловч ерөнхийдөө Эйлэрийн аргаас (тэр үед thetha=1 байсан) илүү багассан байна. Хэрвээ коэффициент тогтмол бол шийдэл (3.8) – аас дахин тодорхойлж болно.

Буржгар хаалтанд байгаа хүчин зүйл нь тогтворжилтон дахь абсолют утга нэгж байдлаас бага байх хэрэгтэй ба үүнийг өсөлтийн хүчин зүйл (amplification factor) гэж нэрлэдэг. Дараах нөхцөлд delta t – ийн ямар ч утганд өсөлтийн хүчин зүйлийг шалгахад хялбархан байх болно.

Үүний чухал зүйл нь хугацааны алхам тогтворжилтоор хязгаарлагдахгүйд орших ба тэгшигэл 3.10 – аас тайвшрах хугацааг 3.4 – д тодорхойлсон шиг дахин олж болно. Бүрэн дифференциал тэгшитгэлийг бодох маш олон арга байх ба бүгд л ил болон далд аргатай байдаг. Зарим нэг ном үзээ.

Gear, C.W. (1971) Numerical Initial Value Problems in Ordinary differential Equations. Prentice Hall

Фортран эмхэтгэл ашиглах учраас фортран хэл дээр төгсгөлөг ялгаварын ил ба далд аргуудыг харьцуулах програмыг дараах байдлаар бичье.

! Computational Hydraulics chapter 2 Section 3

! Water quality in Lake

! Section 3.4 Explicit and implicit

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

          dt=3600 ! [sec] time interval

      nt=500       ! time step

      Q=1.0        ! [m3/sec] Qin, Qout

      V=6.e+5     ! [m3]

      Tr=3.3e+5  ! [sec] = 3.3day

      ci=5.0        ! [mg/l]=[g/m3]

      ds=100.0    ! [g/sec] L discharge

      ce0=(ds+Q*ci)/(Q+V/Tr)    ! eq(2.4)

      T=V/(Q+V/Tr)                    ! eq(2.6)

      write(*,*)"read theta"

      read(*,*)theta

      write(*,'(f9.1,f9.1)') ce0,T/(24.*3600) !; pause

      c0=5.0                               ![mg/l]

      c1=c0                                 ! initial concentration

      c3=c0                                 ! initial concentration

      do n=1,nt                           ! advance of time

      time=n*dt                          ! time

      if(time/(24*3600).gt.3.0) goto 666

! Explicit method

      rh=c1+(dt/T)*(ce0-c1)       ! right hand side of eq.

      c2=rh                                 ! temporal integration

      c1=c2                                 ! next time step

! Mixed method /implicit but theta can give by handle/

      ra=dt*(ce0-c3)-c3*(theta-T)           ! right hand side of eq.3.8

      c4=ra/(T+theta*dt)                         ! temporal integration

      c3=c4                                 ! next time step

! Implicit method

!      ra=dt*(ce0-c3)-c3*(0.2-T)            ! right hand side of eq.3.8

!      c4=ra/(T+0.2*dt)              ! temporal integration

!      c3=c4                               ! next time step

      write(* ,'(4f9.3)') time/(24*3600),c1/ce0,c3/ce0,dt/T! output

      write(10,'(4f9.3)') time/(24*3600),c1/ce0,c3/ce0,dt/T! output

      enddo

      close(10)

  666 stop

      end

     

Figure 3. Тооцоолон бодох аргын үр дүнгийн харьцуулалт