 |
| Төгсгөлөг ялгавар дээр суурилсан даралтын аргаар НС тэгшитгэл болон Конвекци-Тархалтын тэгшигэлийг шийдэх Фортран кодын хэсэг. Урсгалын бусад параметрүүдийг тодорхойлох: Ёроолыг шүргэх хүчдэл, хуйлралт, эсэргүүцлийн коэф CD, Энерги гэх мэт |
Тооцон бодох кодын боловсруулалт
Төгсгөлөг ялгаварын аргын боловсруулалт
Хагшаас ба
давсны тархалтын жишээ болгон давсны шаантагтай голын адгийн хоёр хэмжээст
хүрээг авч үзье. Бодлогийн хүрээг дараах зургаар харуулав.
Figure 4. Бодлогийн
геометр ба хүрээ
Бодлогийн хүрээ нь хязгаарын
нөхцлийн хүрээн дотор байдаг нь загварчлагчийн нэг төрлийн арга юм. Учир нь
ашиглаж буй хязгаарын нөхцөл ихэвчлэн дифференциалчлагдахгүй байхаар өгөгдөх
учир энэ нөхцлийн нөлөөллийг бодлогийн хүрээнээс хол байлгах үүднээс хязгаарын
нөхцлийг гадуур нь авдаг байж болох юм. Ихэвчлэн тусгаарлах хязгаарын нөхцлийг (insulated
boundary condition-Robin b.c) гаралтын хэсэгт хэрэглэх ба бодисын концентраци
болон дулааны тоо хэмжээний градиент хязгаар дээр оршихгүй гэж үзнэ.
С нь ямарч бодисын концентраци болон хэмжигдэхүүний тархалт
байж болох ба дээрх тэгшитгэлийг төгсгөлөг ялгаварын ухрах, давших, төвийн
ялгавараар бичвэл:
Хэрэв cn нь бидний бодлогийн хязгаарын тор (IMAX, JMAX) бол давших болон
төвийн ялгаварын аргыг хэрэглэж болохгүй ухрах ялгавар дээр тунаж үлдэнэ. Иймд
хязгаарын нөхцөл нь хялбар болж
Бодлогийн гол зорилго нь хагшаасны урсгал
болон давсны шаантагын пронт хоёрын харилцан үйлчлэлийг судлахад чиглэж байгаа
юм. Төгсгөлөг ялгаварын аргаар тооцон бодох алгоритмыг авч үзье.
Тооцон бодох алгоритм
Алгоритм гэдэг нь үйлдлийг гүйцэтгэх дэс
дарааллыг тайлбарлах хэрэгсэл бөгөөд тайлбарлах гурван төрлийн арга байх ба
тэдгээрийг блок схемийн (Flow chart), заримдаг кодын (pseudo code), өгүүлбэрийн
(description sentences) гэнэ. Гэхдээ өгүүлбэрийн аргыг алгоритмыг илэрхийлэх
арга гэж үздэггүй, заримдаг кодыг тухайн алгоритмыг боловсруулах хэлний онцлог
дээр тулгуурлаж үзүүлдэг учир өргөн хэрэглэх нь бага байдаг гэж мэт. Тиймээс
заримдаг код нь бодит кодоос ялгаагүй зөвхөн нарийвчилсан үйлдэл командуудыг
хассан илэрхийлэл юм. Блок схемийн аргаар тооцон бодох алгоритмыг зураг 5-аар
үзүүлье. Бодлогийн геометрийг dx=1.0, dx=0.2м-ийн хэмжээтэй тэгш өнцөгт торонд
хуваах ба хэвтээ чиглэлт imax=202 тор, босоо чиглэлд jmax=102 тор байна. Бодлогод
хамаарах фазуудын ялгаатай байдлыг ус, хөрс, давстай ус гэх мэтээр танихын тулд
ялгах параметр ic(i,j)=1,0,2 тус бүрчлэн өгнө. Хэрэв ic параметр 1-тэй тэнцүү
бол түүнийг ус гэж ойлгоно. Ингээд алгоритм дээр суурилж ямар үйлдлийг аль
томъёогоор хэрхэн гүйцэтгэхийг авч үзье.
Анхнынөхцлүүд ба бодлогийн хүрээ
Бодлогийн геометрийг хуваасан торны дагуу
фазыг таних параметр ic – ээр ялгасан файлаар дүрсэлнэ. Геометр хялбар
хэлбэртэй учир програм бичиж үүсгэхэд төвөг байхгүй. Уг файлыг дуудаж ялгах
парамерийг тооцон бодох кодонд уншиж авна.
open(21,file='domain.in')
do 6 j=jm1,2,-1
6 read(21,'(200i1)')
(ic(i,j),i=2,im1)
close(21)
Оролтын хэсэгт хувийн зарцуулгын утгыг q=q+dely*UI
гэж бодох ба оролтын хэсгийн торны индекс нь i=2; j=2,101 хүртэл байх ба dely=dy,
UI нь оролтын хэсгийн хурд 4.5м/с байна. Нийт 100 ширхэг торны хувьд
нийлбэрчилж зарцуулгыг тодорхойлох кодыг бичвэл
q=0.0
i=2
do 7 j=2,jm1
if(ic(i,j).ne.0) q=q+dely*UI
7 continue
Оролтын зарцуулга нь хоёр дахь босоо
нүднүүд дээр тодорхойлогдох ба үүнийг ашиглан i=1 босоо нүдэн дээрх оролтын
хурднуудыг тодорхойлж u(i,j) гэсэн хэмжигдэхүүнд хуулж хадгална. Бусад бодлогийн хүрээн доторхи устай нүднүүд
дахь анхны хурдны утгыг хувийн зарцуулгаас хөөж тодорхойлно.
Figure
5. Тооцон бодох алгоритм /англи ба монгол/
Энд өгсөн хурднууд нь хэвтээ хурд
бөгөөд хөрсний хувьд хэвтээ хурдыг заавал тэг байна гэж өгөх ёстой. Бүх нүдэн
дэх босоо анхны хурдууд нь бүгд тэг байна. Хурдны анхны нөхцлийг өгсөний дараа
нягтын нөхцлийг өгнө. Баруун захад буюу i=imax –аас давсны нягтын үйлчлэл
тэжээгдэх учир энэ баганын нүдэнд i=201-р баганын нүднүүд давстай байх юм бол
j=2,101-ээр нүднүүдэд нягт нь давсны нягт ro=0.01гр/см3, мөн i=imax=202
дугаар нүднүүдэд давс байвал нягтыг ro=0.01гр/см3 гэж өгөөд ялган
таних параметрийг ic(i,j)=2 байсаныг ic(i,j)=1 болгон шилжүүлнэ. Ингэснээр
цаашдын тооцоонд нягт өөр боловч ус хэмээн танигдах болно. Бусад нүднүүдийн
хувьд нягтыг өгөөгүй учир автоматаар хэмээн тооцоонд орно. Энэ нь усны нягт
бөгөөд давстай ус цэвэр уснаас 0.01гр/см3- аар илүү нягттай бодлогод
өгөгдлөө гэж үзэж болно. Шингэний нягтын анхны утгуудыг өгсөний дараа
гидростатикийн даралтыг нүд бүрт бодох болно. Учир нь НС-ийн тэгшитгэлд даралт
үл мэдэгдэгч бөгөөд үүнийг гидростатикийн тэгшитгэлээс олохоос өөр арга
байхгүй. Гидростатикийн даралтыг хамгийн дээд талын нүднээс доош нийлбэрчлэх
замаар бодно. Үүний тулд даралтыг i=2,201 хүртэл, j=jm1=101 дугаартай нүднүүдэд
эхлэн тодорхойлно.
DO 2 I=
2,IMAX
j=jm1
! pressure at
top cell
sum=0.5*dely*ro(i,j)/rowb*(-gy)
P(i,j)=sum
DO 2
J=JM2, 2,-1
! hydrostatic
pressure is calculated by the weight of water over the cell.
sum=sum+(ro(i,j+1)+ro(i,j))/rowb*(-gy)*(0.5*dely+0.5*dely)/2.0
P(i,j)=sum
! pressure in
obstacles = 0.0
if(ic(i,j).eq.0.and.i.ne.imax) p(i,j)=0.0
2 continue
Даралтыг НС дахь хэлбэрээр тодорхойлохдоо дараах ойролцооллыг хийнэ.
Даралтын анхны
утгыг бүх нүдний төвд бодож дууссанаар анхны нөхцлийн дэд програм дуусна.
Бодлогийн хүрээн дотор даралтыг бодсон ба зарим хязгаарын нөхцлийн утга
орхигдоно. Дэд програмын төгсгөлд keeper.out файлыг нээв. Энэ нь хэвлэх бүрт
хэвлэсэн дугаар, хугацааг бичиж хадгалах файл юм.
Хурдны шинэчлэл
Анхны нөхцлөөс
хэвтээ ба босоо хурдны бүрдүүлэгчүүд тус бүр u(i,j); v(i,j) гэсэн хувьсагчаар
хадгалагдсан байгаа. Үүнийг бүхэл нүдний дагуу un(i,j); vn(i,j) гэсэн
хувьсагчид хуул дараагийн алхамд НС болон конвекци тархалтын тэгшитгэлд
хэрэглэхэд бэлдэнэ.
Тоон шийдлийн тогтворшилт, хугацааны алхамын сонголт
Тогтмол хугацааны
алхамтай бодолт нь шийдийн тогтворгүй байдалд хөтөлдөг. Шийдийн тогтворшилт,
хугацааны алхамд шинжилгээ хийх Вон Неуманн, CFL нөхцлүүд зэрэг байдаг. Эдгээр
шинжилгээний аргууд нь хурдны градиент болон орон зайн алхамтай хугацааны алхам
хэрхэн хамаарч байгааг тодруулдаг. Нэгэнт орон зайн алхамыг сонгосон учир түүн
дээрх тэгшитгэлийн шийдийн нийлэмж, дөхөлтийг хугацааны алхамаар зохицуулах нь
илүү хялбар. Коурантын нөхцөл ёсоор хугацааны алхам нь дараах нөхцлийг хангаж
байх ёстой.
Дээрх илэрхийллийг
хувирган нягтыг харьцангуй тогтмол гэж үзэн нэгтгэн бичиж хугацааны алхамыг
сонгох боломжийг хайвал:
Хоёр хэмжээст
бодлогод тохируулан тогтворшилтын нөхцлийг бичвэл:
Конвекци-тархалтын
нөлөө орох учир хурдын хамгийн их утгаар хугацааны алхамыг үнэлэх нь оновчтой
ба мөн нэмэлт хязгаар тавьж хугацааны алхам ямагт 0.1сек – ээс бага байхаар
хязгаар тавих нь нэг нөхцөл юм. Ингээд тооцооны бодит хугацаа болон үр дүнг
хэвлэх хугацааг дараах байдлаар тодорхойлно.
Үр дүнг tprt =0.5сек давтамжтайгаар
хэвлэх (np1=0.5, np2=1, np3=1.5, np4=2
гэх мэт
хугацаанд үр дүнгүүд хэвлэгдэх ба энэ нь бодит хугацааг илтгэнэ) ба нэг давтамж тодор сонгогдсон хугацааны
алхамаас хамаарч хэдэнч цикл тооцоо хийгдэж болно. Иймд нийт циклийн тоог
хүрэлцээтэй байхаар өгөх ба жишээнд nfinl=1’000’000 гэж өгсөн байна.
Ингээд бодлогийн гол
цикл эхлэнэ.
Үндсэн гогцоо
Үндсэн циклийн эхэнд
тооцооны бодит хугацаа үр дүн хэвлэх хугацаанд хүрч ирсэн бол үр дүнг хэвлэнэ.
Үр дүнг хэвлэх байрлал алгоримтын хаана байрлах нь заримдаа чухал байдаг. Хэвлэх
дэд програмын тухай сүүлд тайлбарлана.
Нягтын үйлчлэл ба хагшаасны концентрацийн зөөгдөлт-тархалт
Энэ дэд програмд
бодлогийн геометрийн зүүн талаас хагшаасны концентраци тэгшитгэл 22- ийн дагуу,
баруун талаас нягтын үйлчлэл тэгшитгэл 21-ийн дагуу тус бүр тархах болно. Тэгшитгэлүүдийг
төгсгөлөг ялгаварын аргаар тухайлж бичвэл:
Үүнтэй адилаар хагшаасны
концентрацийн тэгшитгэлийг бичвэл:
Дээрх
тэгшитгэлүүдийн зөөгдөлтийн гишүүнийг эсрэг ялгаварын схемээр тухайлан бичье. Доорхи тэгшитгэлд frx, fry гэх мэт нь зөвхөн
ялгах таних тэмдэглэгээ юм.
Зөөх хурдны
тэмдэгээс хамаарч эсрэг ялгаварын схемийг тэгшитгэл 40, 41-ийн дагуу бодно. Нягт
ба концентрацийн тархалтын гишүүдийг дараах байдлаар тухайлна. Ингэхдээ
хоёрдугаар эрэмбийн дифференциал учир төвийн ялгаварын схемийг ашиглана.
Хязгаар дээр
ялгаварын гурван зангилаа ашиглагдах боломжгүй учир хэвтээ ба босоо чиглэлд дахь төвийн ялгаварыг
задалж i ба j зангилааны хоёр талд бодно. Тухайлбал:
Босоо тэнхлэгийн
хувьд ижил замаар бичигдэх ба эцсийн бүлэгт дээрх тус бүр бодогдсон гишүүдийг
тэгшитгэл 38, 39-д орлуулсанаар нягтын үйлчлэл ба хагшаасны концентраци
бодогдож дуусна.
Навьер-Стокесийн тэгшигэлийг шийдэх
Хоёр хэмжээстэд
тэгшитгэл 26-г задлан бичиж төгсгөлөг ялгавараар ойролцоолъё.
Дээрхи тэгшитгэлийг
шийдэх нь энэ дэд програмын зорилго байх болно. Шингэний нягтыг өмнөх
конвекци-тархалтын тэгшитгэлээр бодсон үр дүнгээс, даралтыг анхны нөхцөл дахь
гидростатикийн даралтаар авна. Дагаварт тэмдэглэгээтэй байгаа адвекцийн
илэрхийлэлийг босоо ба хэвтээ орон зайд хоёр хувааж тодорхойлох ба ойролцоололд
нэгдүгээр эрэмбийн эсрэг болон төвийн ялгаварын хос схемийг хэрэглэнэ. Төвийн
ялгаварын схем нь шийдийн нарийвчлал сайн байх боловч тогтворжилт багатай бол
нэгдүгээр эрэмбийн эсрэг схем нь тогтворшилтыг хангаж өгнө. Зунгааралтын
гишүүнийг жишээ болгон хэвтээ тэнхлэгийн дагуу ойролцоолбол:
НС- ийн тэгшитгэлээр
урсгалын шинэ хурдны бүрдүүлэгчүүд бодорхойлогдоно.
Шийдлийн нарийвчлалыг шалгах гогцоо
НС- ийн тэгшитгэл,
конвекци тархалтын тэгшитгэл зэргийг тооцон бодож хурд болон бусад физик
хэмжигдэхүүнийг тодорхойлсон ба уг бодлогийн нөхцөл нь үл шахагдах шингэний
хөдөлгөөн учир олсон хурд нь нөхцөл хангаж байгаа эсэхийг тасралгүйн
тэгшитгэлээр шалгана. Босоо ба хэвтээ хурдны градиентуудын нийлбэр тэг байх
ёстой боловч EPSI=0.01-ээс бага байх шаардлагыг тавья. Хэрэв нөхцөл
хангагдахгүй бол физик хэмжигдэхүүнд алдааг хуваарилаж нөхцөл хангагдах хүртэл
итераци хийнэ. Итерацийн хамгийн их хэмжээ 500 удаа байна. Бодолтын дөхөлтийг
шалгахын өмнө хязгаарын нөхцлийг өгнө.
Итераци дахь хязгаарын нөхцөл
Оролт ба гаралтын
хэсэгт бодлогийн хүрээн дотор тодорхойлогдсон нягт, концентраци ба хурдны
бүрдүүлэгчийг хязгаарын хүрээний утгуудруу хуулна. Тооцооны онлогоос хамаарч
ханын нөхцлийг гулсах болон үл гулсах байдлаар өгч болно. Ханын материалын
үрэлтийг эс тооцох бол гулсах нөхцлийг өгч болох ба энэ нөхцөлд хана орчидм
хурдны градиент үүсэхгүй. Үрэлтийг тооцох бол үл гулсах нөхцлийг өгнө. Жишээ
болгон далайн ёроолд гулсах ба үл гулсах нөхцлийг өгье.
if(ic(i,j).eq.0.and.ic(i+1,j).eq.0.and.ic(i,j+1).eq.1) then
if(ibc.eq.1) then
! ibc=1 slip boundary
condition on wall
U(I,J)= U(I,J+1)
else
! ibc ne.1 no slip boundary
condition on wall
U(I,J)=-U(I,J+1)
endif
endif
Гулсах ба үл гулсах
нөхцлийн тодорхойлгоч параметр нь ibc бөгөөд үүнийг оролтоор хянах болно. Авч
үзэж буй зангилааны дараачийн зангилаа нь хөрс, дээд зангилаа нь ус байх
тохиолдолд уг зангилаан дээрх хурдыг
1.
Үл гулсах нөхцөлд uij=-ui,j+1 гэж
өгнө.
2.
Гулсах нөхцөлд uij=ui,j+1 гэж өгнө.
Тоон шийдийн дөхөлт шалгах
Тоон шийдлийн
дөхөлтийг зангилаа бүрд шинжилнэ. Зангилаан дахь градиентын хазайлтыг ухрах
ялгавар ашиглан тэгшитгэл 2-оос дараах байдлаар тодорхойлно.
Зангилаа
бүрийн утгын хазайлтыг бусад зангилаатай харьцуулж аль их байгааг нь тараах
алгоримт үйлчлэх ёстой. Тараах хэмжээг дараах байдлаар тодорхйолно.
Үл шахагдах нөхцөл
хангаж шийдэл нийцэлтэй болсон бол алгоритмын дараагийн дэд програмд шилжиж
шинэчлэн бодогдсон хурдыг дараачийн тооцооны алхамд шилжүүлнэ. Үүний дараа
дахин тооцооны тогтворшилтын нөхцлийг оруулж шинэ хугацааны алхамыг бодож
тооцооны 2-р циклдээ орно. Цикл эхлэхэд мөн л хэвлэх хугацааг шалгаж логик
биелэвэл үр дүнг хэвлэх дэд програмыг дуудна.
Шингэний урсгалын бусад хэмжигдэхүүнийг тооцоолох
Хагшаас ба давсны
түрлэгийн харилцан үйлчлэлийн улмаас хуйлралт ажигдагдах ба үүнийг хуйлралтын
хялбаршуулсан тэгшитгэлээр торорхойлж болно.
 |
| Навиер-Стокесийн тэгшитгэлийг шийдэх хэсэг. Код хэт урт учраас алгоритмлж тус бүр дэд програмд ямар үйлдэл хийгдэх, ТЯ-аргаар хэрхэн дифферецниал тэгшитгэлийг ойролцоохлох гэх мэтийг үзүүлэв. Үр дүнг 10 сард оруулна. |
No comments:
Post a Comment