Бүлэг 7

Тархалт буюу диффузи /diffusion/

7.1 Хэвтээ үе дахь газрын доорхи усны урсгал

Хөрсий сүвэрхэг хэвтээ үе дахь газрын доорхи усны урсгалыг авч үзье (зураг 7.1). Хөрсний босоо багана дахь массын балансын тэгшитгэл нь голын массын балансын тэгшитгэлтэй адил байна.

Энд a=уст үеийн зузаан, u=үеийн зузаанд дундажилсан урсгалын хурд, w=нэгж талбайд орох хур тунадасны эзлэхүүн, n=хөрсний сүвэргэшил

Борооны ус нь газрын доорхи усны түвшинг нэмэгдүүлнэ. Энд хөрсний дээр байх усаар ханаагүй хөрсний үеийг авч үзээгүй. Хэрэв авч үзвэл тухайн үе усаар ханасаны дараа газрын доорхи усны түвшин нэмэгдэх болно. Газрын доорхи усны урсгалын динамик тэгшитгэл болох Дюпюьт-ийн тэгшитгэлийг авч үзье.

Үүнд: h=хэвтээ хавтгайтай харьцуулсан газрын доорхи усны түвшин, k=хөрсний нэвчүүлэх чанарыг илэрхийсэн тогтмол (энэ нь хөрсний шинж чанарыг харуулах ба нэгж нь м/с байна.).

Дээрх хоёр тэгшитгэлийг нэгтгэвэл

Дээрх тэгшитгэл нь энэ бүлгийн гол тэгшитгэл юм.

Олон тохиолдолд, гүний усны үеийн зузаан a их хэмжээгээр өөрчлөгддөгүй учир D=ka нэгтгэл нь ойролцоогоор тогтмол байна (түүгээрч барахгүй нэвчилт k нь ихэвчлэн сайн тодорхойлогдоогүй байдаг). Ингээд стандарт тархалтын тэгшитгэлийг гаргаж авбал

Энд D нь тархалтын коэффициент ба нэгж нь м²/с байна. Энэ төрлийн тэгшитгэл физик болон математикт жишээлбэл дулааны дамжилт гэх мэтээр нилээн тохиолдоно. Энэ бол x –ийн хувьд хоёрдугаар эрэмбийн, t – ийн хувьд нэгдүгээр эрэмбийн дифференциал тэгшитгэл юм. Энэ шинж чанараараа л өмнөх бүлэгт дурьдагдаж байсан энгийн долгионы тэгшитгэлээс үндсэндээ ялгаатай юм.

Голын хажуу талд гүний усны хэвтээ давхрагын жишээгээр тархалтын үзэгдлийг тухай санааг авъя. Бүх зүйл тэнцвэрийн нөхцөлд байна гэж үзье (хур тунадасгүй, газрын доорхи усны түвшин хэвтээ байдалтай). Ямар нэгэн нөхцөлийн улмаар голын усны түвшин доошилж (газрын доорхи усны хариу үйлдэлтэй харьцуулахад маш түргэн) тогтмол болно. Энэ үед газрын доорхи ус голруу урсаж зураг 7.2 – т үзүүлсэнчлэн гадаргын хэвгий (тэгшитгэл 7.2) үүснэ.

Энд erfc гэдэг бол стандарт математикийн функ (төгсжилөртийн алдааны функц) юм. Энэ нь дараах байдлаар тодорхойлогдоно.

Энэ функц z нь хязгааргүйрүү тэмүүлж байхад 0 болох ба z=0 болоход нэг гэсэн утгыг өгнө. 

Үүний эгшин зуурын хугацаан дахь шийдлийг зураг 7.3 – т үзүүлэв. Тэгшитгэл 7.5 болон зураг 7.3 – аас тархалтын ердийн шинж чанаруудыг дүгнэж болохоор байна.

        i.            Ямар нэгэн t хугацааны дараа, хэдий энэ өчүүхэн бага байсанч, x=0 үед дахь хуваарилалтын нөлөө нь хаана ч тэр мэдрэгдэнэ (сулхан ч мэдрэгдэж болно). Хуваарилалт нь хязгааргүй хурдаар өөрөөрөө өргөжнө. Газрын доорхи усны бууралт нь x/t^1/2=const үед ижил байна. Тиймээс хуваарилалтын өргөсөлт нь хугацааны хувьд шугаман бус гэвч t^1/2 шиг л байна. Өргөсөлт нь мөн D^1/2 – той шууд хамааралтай байна. Нөлөөллийн ач холбогдол нь x=(4Dt)^1/2 – ийн утганд мэдрэгдэнэ.

    iii.            Голоос L зайнд өөрчлөлтийг мэдэх боловч газрын доорхи усны түвшний бодит өөрчлөлт нь L²/4D –ийн эрэмбийн хугацааны дараа ажиглагдах болно.

Анхын болон хязгаарын нөхцөл нь энгийн долгионы тэгшитгэлийг бодвол жаахан хэцүү гэвч дараах итгэл төрүүлэхүйц дүрмүүдийн дагуу тодорхойлж болох юм (тухайн дифференциал тэгшитгэлийн ямарч сайн номонд энэ тухай дэлгэрэнгүй байгаа, жишээ нь Garabedian, 1967). Тархалтын тэгшитгэл нь хугацааны хувьд нэгдүгээр эрэмбийн, тэгэхээр нэг анхны нөхцөл хэрэгтэй, ихэвчлэн t=0 үе дахь газрын доорхи усны түвшин оновчтой. Бидэнд u- ийн тодорхой утга байхгүй байж болно, үүнийг тэгшитгэл 7.2 – оос h – ийг ашиглан гаргаж болно.

Зайн дахь хоёрдугаар эрэмбийн тэгшитгэлд хоёр хязгаарын нөхцөл хэрэгтэй. Аль аль талд тус бүр нэг.

-          Голын захад (x=0) нийлэх газрын доорхи усны түвшин нь голын усны түвшинтэй тэнцүү гэж авъя.

-          Хязгааргүйд орших (Голоос хангалттай хол орших өөрөөр хэлбэл голоос хамгийн бага зайнд) газрын доорхи усны түвшин тогтмол гэж авъя.

7.2 Ил төгсгөлөг ялгаварын арга

Хэрэв тархалтын тэгшитгэлийг тоон аргаар шийдэх гэж байгаа бол эхлээд хоёрдугаар эрэмбийн гаргалгааг ойролцоогоор олох хэрэгтэй. Энэ нь нэгдүгээр эрэмбийн гаргалгааг давтах процессоор хийгдэнэ.

Дээрх бодолт нь төвийн төгсгөлөг ялгаварын ойролцоолол юм. Үүнийг урагш хугацааны ялгавартай нэгтгэвэл урагш хугацааны, төвийн зайн эвсэл FTCS аргыг гаргах болно.

Энэ тэгшитгэл нь өгөгдсөн анхны нөхцлөөс эхлэх болно. Хэрэв, хязгаарын нөхцөлд h хязгаарын нөхцлөөс мэдэгдэж байвал, ямарч онцгой гаргалгаа шаардлагагүй. Гэсэн хэдий ч хэрэв dh/dx өгөгдсөн бол салангад байх шаардлагатай.

Одоо тогтворжилтгүй байдлын талаар санаа зовох ёстой учир тогтворжилтыг шалгаж болно. Нэмэгдэлтийн хүчин зүйл нв хялбархан дараах байдлаар олдоно.

Энд шинэ хэмжээсгүй параметр бол lamda=2Ddelta t/delta x² болох ба үүнийг тархалтын параметр гэж нэрлэнэ. Нэмэгдэлтийн хүчин зүйл нь бодит учир ямарч psi – ийн абсолют утганд нэгжээс хэтрэхгүй болохыг харуулахад хялбархан.

Энд дахин хугацааны алхамын “торны хэмжээн квадраттай” холбогдсон хязгаарлалтын тухай асуудал бий болно. Хэрэв нарийвчлалтай байх үүднээс жижиг хэмжээтэй тор ашиглах гэж байгаа бол муу байх болно.

Жишээнд, дараах утгуудыг авъя. D=10^-3м/с, delta x=10м, delta t=15цаг.

Усны түвшний буурал delta h – ийн утга нь тархалтын тэгшитгэл шугаман учир ач холбогдолгүй. Зураг 7.4а нь 10, 20, 30 хугацааны алхамын дараах шийдлийг харуулж байна. Зургийг харвал, энэ тохиолдол нь тогтвортой биш ба энэ нь lamda=1.08 утганд батлагдаж байна. Хэрэв хугацааны алхамыг хоёр дахин багасгавал зураг 7.4b харуулсан шийдэл гарах ба энэ нь тогтвортой байна.

Мэдрэхүй дэхь тархалтын параметр lamda нь Коурантын тоотой харьцуулагдаж болохоор байна. Энэ нь L²/4D эрэмбэтэй хугацаанл авсан L уртын дагуух тархалтын хуваарилалтын нөлөөг дурьдаж байна. Delta x гэсэн нэг торны интервалд энэ нь delta x²/4D байна. Тогтворшилтын нөхцөл (7.9) нь хугацааны алхам “тархалтын хугацаа”наас хэтрэхгүй байх ёстой гэдгийг харуулж байна.

7.3 Далд төгсгөлөг ялгаварын арга

Тэгшитэл 7.9 дахь тогтворжилтын хязгаарлалтыг далд арга ашиглан давж гарах зам байдаг. Санаа нь гэдэл аль хэдийн 3.4 хэсэгт өгөгдсөн байгаа. Гол зарчим нь хугацаа t –д биш ямар нэгэн t+thetha delta t хугацаанд (0.ge.thetha.le.1 байна.) ялгаварын тэгшитгэлээр ойролцоолох явдал юм. Тэгвэл

Гэвч орон зайн гаргалгаанд бид хоёр хугацааны түвшний дунджыг авах хэрэгтэй ба салангад торонд ямар ч дундажилсан утга байхгүй.

Дээрхийг бүгдийг нь тархалтын тэгшитгэлд орлуулбал үр дүнд нь Кранк-Ничолсоны арга гарч ирнэ.

Үүний хэцүү нь юу вэ гэвэл n+1 хугацааны түвшин дахь гурван үл мэдэгдэх утга агуулсан явдал юм. Тэгэхээр энэ өөрөөрөө шийдэгдэх боломжгүй гэсэн үг. Гэсэн хэдийч тэгшитгэлийг тоолбол, дотоод цэг бүрт тэгшитгэл 7.12 нэг, хоёр хязгаарын цэг бүрт нэг хязгаарын нөхцөл байгааг харах болно. Тиймээс, бүх үл мэдэгдэгч нь системээс шийдэгдэх боломжтой байна. Энэ нь тэгшитгэлээс нэг ил утга олох боломжгүй гэвч тэдгээрийг далд байдлаар олох боломжтой байна. Энэ л энэ төрлийн аргын нэрний тайлбар юм. Ижил санаа энгийн долгионы тэгшитгэлийн шийдэлд хэрэглэгдсэн болно. Тэгшитгэлийн системийн шийдвэрлэх дээд зэргийн аргыг 7.4 хэсэгт тайлбарлах болно. Далд шийдлийн жишээ нь зураг 7.5 – д өгөгдсөн байна. Энэ нь зураг 7.4 – тэй ижил нөхцөлд холбогдоно. Энд thetha=1/2 ба delta t=28цаг эсвэл 140цаг нь lamda=2 эсвэл 10 утгатай тохирч байна. Альч шийдэл төгс тогтортой байх ба нарийчлал нь ижил биш байна.

Энэ арга нь lamda>1 утганд үнэхээрийн тогтвортой байгаа нь нэмэгдэлтийн хүчин зүйлийг хэрэглэсэнээр харагдана.

Үүнийг ашиглан thetha – ийн доорхи утгатай байхад lamda- ийн ямарч утганд энэ аргыг тогтвортой болохыг энгийнээр харуулж байна.

Та тархалтын параметрийн ач холбогдол дээр 7.2 хэсэгт тайлбарласантай санал нэг байсан ч гайхаж магадгүй. Энэ тохиолдолд j цэг дээрх үл мэдэгдэх утга өөрийн хөршдөө нөлөөлөх ба эдгээр нь эргээд өөрсдөдөө улмаар хязгаарын цэгт нөлөөлнө. Бас өмнөх хугацааны түвшний бүх утгууд өөрсдөдөө нөлөөлж болно. Тиймээс хуваарилалт нь нэг хугацааны алхамын турш бүхэл бүсэд нөлөөлөлтэй байж болно гэдэг дифференциал тэгшитгэлийн шиж чанарыг зөвшөөрөх хэрэгтэй. Энэ нь хугацааны алхамыг ямарч утгад үнэн байна.

7.4 Томасын алгоритм

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

Гол диагнал болон түүний хажуу хөршүүдээс бусад матрицын гишүүдийн ихэнх нь 0 байна. Үүнийг гурван диагональт матриц гэнэ.

Энэ төрлийн системийн шийдэл нь Томас эсвэл давхар шүүх алгоритмаар маш үр дүнтэйгээр шийдэгддэг. Давхар шүүх алгоритм гэдэг нь стандарт Гауссын зайлуулах аргатай ижил.

Хэрэв дараах хэлбэрт ямар нэгэн тэгшитгэл бичигдсэн байвал

а_o=c_J=0 үед дараах нөхцөл бүрдэнэ.

Энэ нь дараах зүйлийг гаргана.

Үүнийг j-1 дахь тэгшитгэл 7.16 –г тэгшитгэл 7.15 – д орлуулсанаар нотолж болно. Тэгшитгэл 7.17 болон 7.18 –ын эхнийх нь j=1,2,…….,J – ийн хувьд шийдэгдэнэ (урагш шүүх арга). Тэгээд h_J – г хязгаарын нөхцлөөс тодорхойлно. Төгсгөлд нь арагш шүүх аргаар бусад үл мэдэгдэх гишүүдийг тэгшитгэл 7.16 аас j=J-1,J-2,…….,1 гэх мэтээр тодорхойлж болно. Процесс нь компьютерийн хугацаа болон багтаамжид үр ашигтай байдаг.

Гурван диагональт матрицийн систем

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

Үүнийг n хэмжээст гурван диагональт матриц хэлбэрээр бичвэл

Энэ матрицийг бодох алгоримтыг бичвэл (Гауссын зайлуулах арга)

do k=2, n

m=a_k/b_k-1

b_k=b_k-m*c_k-1

d_k=d_k=m*d_k-1

end do

x_n=d_n/b_n

do k=n-1, 1

x_k=(d_k-c_k*x_k-1)/b_k

Манай дээрх 7.4 дэхь жишээний хувьд бичхээс залхуу хүрсэн учир зургийг нь оруулчий.

7.5 Хэрэглээ

Хөдөө аж ахуйн шүүрүүлийн талбайд хоорондоо L зайтай хэвтээ ил шүүрүүл байна гэж үзье. Т хугацааны турш q=w/n эрчимтэй бороо орж гэж үзье. Ингэхэд газрын доорхи усны түвшин хэр өргөгдөж, ямар хугацааны дараа дахин хэвийн байдалд орох бол?

Хөрсийг нэгэн төрлийн гэж үзье. Мөн сувгийн хоорондох нөхцөл нь тэгш хэмт учир L/2 өргөнтэй талбайг л авч үзэхэд хангалттай. Тэгш хэмийн шугаманд хязгаарын нөхцөл dh/dx=0 байна. Шүүрүүлэх сувагт, усны түвшин тогтмол үлдэнэ гэж үзье (h=0). Үүнтэй адил усны түвшинг анхны нөхцөл болгон авъя.

L=100м гэвэл хагас бүсэд delta x=2.5м – ийн интервалаар 20 тор авъя. Тархалтын коэффициентыг өмнөх жишээтэй ижил байдлаар D=10^-3 м²/с гэж авъя. Бороо орох хугацааг T=29хоног гэвэл хугацааны алхам delta t=0.5өдөр нь хангалттай жижиг мэт санагдана. Тэгээд lamda=15, гэвэл далд арга хэрэгтэй болно. thetha=1/2 утгатай Кранк-Ничолсоны аргыг хэрэглэж зураг 7.6 харуулсан үр дүнг гаргаж авав. Энд усны түвшиний нэмэгдэлтийн нэгжгүй байдлаар харуулав. Бороо орж дуусах хугацаанаас T хугацааны дараа гүний усны түвшин буцаж хэвийн байдалд орохоор байна.

Тархалтын тухай энэ бүлэг үүгээр өндөрлөв. Тооцон бодох жишээгээр арай өөр жишээ авсан болно. Авч үзэж буй ..м газарт усны түвшин гадаргаас доош 10м байх ба хэсэгхэн хэсэгт 9м байна гэж үзээд мөн л тэр хэсэгхэн газарт бороо орж байвал хугацаанаас хамаарч ус гүнрүү хэрхэн нэвчих вэ гэдгийг тооцож үзье.

1.      Ил арга /Explicit method/ Фортран хэл

Computational Hydraulics Sec.7 Diffusion

! computational hydraulics chapter 7 pp.37

! diffusion equation

! equation (7.7)

      parameter(g=9.81,pi=3.141592 )

      parameter(nt=2000,dt=1.0) ! time

      parameter(nx=100,dx=2.0)  ! space

      parameter(difc=1.2,ramda=2.0*difc*dt/dx**2) ! pp.37, I.16

      parameter(rain=0.00001) ! m/s=w/n in the text book

      dimension h1(nx),h2(nx)    ! definition p.34

       do i=1,nx ! initial water level

                      h1(i)=-10.0

      if(i.gt.nx/2-5.and.i.lt.nx/2+5) h1(i)=-9.0

      enddo

      open(10,file='h.dat') ! for data storage

      do n=1,nt ! temporal change

      time=dt*n ! time

      do i=2,nx-1

      if(i.gt.nx/2.and.i.lt.nx/2+20) then

       rh=rain

      else

       rh=0.0

      endif

      rh=rh+ramda/2.*(h1(i+1)-2.0*h1(i)+h1(i-1))+rain

      h2(i)=rh*dt+h1(i)

      enddo ! i

      h2(1)=h2(2)

      h2(nx)=h2(nx-1)

      h1=h2

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

      write(10,'(500e15.4)') (h1(i),i=1,nx)

      endif

      enddo ! time

      close(10)

      stop

      end

Хэрэв усны түвшин бүхэл хэсэгт 10м байсан бол тархалт нь

Computational Hydraulics Sec.7 Diffusion

! computational hydraulics chapter 7 pp.37

! diffusion equation

! equation (7.7)

      parameter(g=9.81,pi=3.141592 )

      parameter(nt=2000,dt=1.0) ! time

      parameter(nx=100,dx=2.0)  ! space

      parameter(difc=1.2,ramda=2.0*difc*dt/dx**2) ! pp.37, I.16

      parameter(rain=0.00001) ! m/s=w/n in the text book

      dimension h1(nx),h2(nx)    ! definition p.34

       do i=1,nx ! initial water level

                      h1(i)=-10.0

      ! if(i.gt.nx/2-5.and.i.lt.nx/2+5) h1(i)=-9.0

      enddo

      open(10,file='h1.dat') ! for data storage

      do n=1,nt ! temporal change

      time=dt*n ! time

      do i=2,nx-1

      if(i.gt.nx/2.and.i.lt.nx/2+20) then

       rh=rain

      else

       rh=0.0

      endif

      rh=rh+ramda/2.*(h1(i+1)-2.0*h1(i)+h1(i-1))+rain

      h2(i)=rh*dt+h1(i)

      enddo ! i

      h2(1)=h2(2)

      h2(nx)=h2(nx-1)

      h1=h2

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

      write(10,'(500e15.4)') (h1(i),i=1,nx)

      endif

      enddo ! time

      close(10)

      stop

      end

2.      Далд арга /Implicit with Thomas algorithm/

Computational Hydraulics Sec.7 Diffusion

! computational hydraulics chapter 7 pp.37

! diffusion equation

! equation (7.7)

      parameter(g=9.81,pi=3.141592 )

      parameter(nt=2000,dt=1.0) ! time

      parameter(nx=100,dx=2.0)  ! space

      parameter(difc=1.2,ramda=2.0*difc*dt/dx**2) ! pp.37, I.16

      parameter(rain=0.00001) ! m/s=w/n in the text book

      dimension h1(nx),h2(nx),a(nx),b(nx),c(nx),d(nx)    ! definition p.34

           do i=1,nx ! initial water level

           h1(i)=-10.0

      !if(i.gt.nx/2-5.and.i.lt.nx/2+5) h1(i)=-9.0 /– Энд тайлбарын тэмдэг ! байхгүй үед эхний үр дүнг авна. Зураг №/

      enddo

      write(*,*) "input thetha"

      read(*,*)thetha

         open(10,file='h3.dat') ! for data storage

      do n=1,nt ! temporal change

      time=dt*n ! time

       do i=1,nx

      if(i.gt.nx/2.and.i.lt.nx/2+20) then

       rh=rain

      else

       rh=0.0

      endif

      rh=rh+ramda/2.*(1-thetha)*(h1(i+1)-2.0*h1(i)+h1(i-1))+rain

      d(i)=rh*dt+h1(i)

       enddo ! i

      d(1)=h1(2)

      d(nx)=h1(nx-1)

      write(*,*) d(1)

     call trid(a,b,c,d,nx)

      h1=d

      h1(1)=h1(2)

      h1(nx)=h1(nx-1)

      !h1=h2

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

      write(10,'(500e15.4)') (h1(i),i=1,nx)

      endif

      enddo ! time

      close(10)

      stop

      end

! Thomas algorithm

      subroutine trid(a,b,c,d,nx)

      dimension a(*),b(*),c(*),d(*)

       do i=1,nx ! matrices

      a(i)=-0.5*thetha*ramda

      if(i.eq.1.and.i.eq.nx) a(i)=0

      b(i)=1+ramda*thetha

      if(i.eq.1.and.i.eq.nx) b(i)=1

      c(i)=-0.5*thetha*ramda

      if(i.eq.1.and.i.eq.nx) c(i)=0

       enddo

      do i=1,nx-1

      m=a(i+1)/b(i)

      b(i+1)=b(i+1)-m*c(i)

      d(i+1)=d(i+1)-m*d(i)

       enddo

      d(nx)=d(nx)/b(nx)

       do i=nx-1,1,-1

      d(i)=(d(i)-c(i)*d(i+1))/b(i)

       end do

      end