Шийдлийн алгоритм дахь Шингэний эзлэхүүний арга (Volume of Fluid Method in SOLA)

Шийдлийн алгоритм дахь (SOLA)

Шингэний эзлэхүүний арга (VOF)

Шингэний чөлөөт гадаргууг ойролцоолох нь тооцон бодох гидравликийн нэгэн чухал асуудал байсаар ирсэн. Ойрын хэдэн жил байх ч болно. Чөлөөт гадаргууг тодорхойлох батлагдсан аргуудыг доорхи зургаас хараарай. Хамгийн анхны арга нь Маркер болон нүдний арга юм. Энэ аргад Лагранжийн координатаар явах маркар буюу цэгнүүд нь шингэний чөлөөт гадаргууг тодорхойлоход хэрэглэгддэг. Үүний дараа гарч ирсэн одооны ихэнх арилжааны програмуудын /Ansys, Solidworks, Xflow, 3Dflow гэх мэт/ ашигладаг арга бол шингэний эзлэхүүний арга /Volume of fluid/ юм. Энэ удаад хамгийн анх боловсруулагдсан шингэний эзлэхүүний аргын тухай авч үзье. Энд тооцооны төгсгөлөг нүдэнд байх шингэнийг эзлэхүүний хувиар илэрхийлж түүнийгээ шингэний фракц гэж авч үзнэ. Шингэний эзлэхүүний фракцийн хугацаан дахь өөрчлөлтийг энгийн адвекцийн тэгшитгэлээр шийднэ.

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

Шингэний чөлөөт гадаргууг ойролцоолох тооцонгийн аргууд.

Шийдлийн алгоритмд НС-ын тэгшитгэлийг ойролцоолох ба

Х чиглэлд

У чиглэлд

гэж бичигдэнэ. Үүний дараа гарч ирсэн хурд болон даралтыг урсгал тасралтгүйн тэгшитгэлийг хангаж байгаа эсэхээр шалгана.

Урсгал тасралтгүйн тэгшитгэлийн ерөнхий хэлбэр

ба хязгаарлагдмал шахагдалтыг тооцвол

болно. Үүнд c нь шингэн дахь дууны хурд, кси нь координатын параметр юм. Кси –г нэг гэж авч үзвэл координатын систем цилиндр болох ба координатууд (x,y)-н оронд (r,z) бичигдэх болно. Харин кси тэг бол тэгш өнцөгт координатын схемд бичигдэх ба 3-5 тэгшитгэлүүдийн сүүлийн гишүүд үгүй болно.

Ерөнхийдөө ША-д шингэний эзлэхүүний аргыг хэрэглэхэд дараах үндсэн гурван алхамыг гүйцэтгэх хэрэгтэй.

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

Дээрх алхамууд нь тодорхой хугацааны тохирол хүртэл хийгдэх ба мэдээж хугацааны алхам бүрт тохирсон хязгаарын нөхцлүүдийг бодлогын хүрээ, чөлөөт гадаргуу дээр өгөх хэрэгтэй.

Моментийн тэгшитгэлийг ойролцоолох нь

Дараах Q_i,jⁿ гэсэн тэмдэглэгээ нь Q(x,y,t) –ийн утгыг илэрхийлэх ба энэ нь хугацаа ndt үед, х чиглэл дахь ith нүд, у чиглэл дахь jth нүдний утга юм. Хагас бүхэл индекс нь нүдний захуудыг илэрхийлнэ. Ингээд тэгшитгэл 2 ба 3-ыг төгсгөлөг ялгавараар ойролцоолсоныг бичвэл:

болох ба үүнд

байна. Эдгээрт нэмэгдсэн шингэний эзлэхүүний утгууд нь хийн нүдэн дахь даралтын утгыг засварлана. Адвекц болон зунгааралтын хурдатгалууд нь FUX, FYX, VISX гэх мэтээр тэмдэглэгдсэн ба тухайлбал FUX гэдэг нь х чиглэл дахь хэвтээ хурд u –ийн адвекцийн цүлхэлт юм. эдгээрийн утгууд нь бүгд хуучин хугацааны алхамын утгуудыг ашиглан олдоно. Яагаад гэвэл шинэ хугацааны алхам дахь даралтууд мэдэгдэхгүй ба тэгшитгэл 6 нь uⁿ⁺¹,vⁿ⁺¹ утгуудыг шууд олоход хэрэглэгдэж чадахгүй тул дараах маягаар урсгал тасралтгүйн тэгшитгэлтэй хавсрана. Эхний алхамд p_n+1  -ийн оронд p_n утга хэрэглэгдэж анхны шинэ хурднуудыг гаргаж авна. Тэгшитгэл 6-ийн адвекц болон зунгааралтын гишүүд нь ач холбогдол багатай боловч тооцонгийн тогтворшилтыг хангахаар бэлтгэгдсэн байна.

Зураг 2. Хэмжигдэхүүнүүдийн байрлал тэгш өнцөгт төгсгөлөг торон дээр.

Урсгал тасралтгүйн тэгшитгэлийг ойролцоолох нь

Тэгшитгэл 6-аас тодорхойлогдсон хурднууд нь урсгал тасралтгүйн тэгшитгэл 5-г хангаж байх ёстой. Энэ тэгшитгэлийг хангахын тулд даралт нь шингэнээр эзлэгдсэн тооцооны нүд болгонд тохируулагдах ёстой. Тэгшитгэл 5-д төгсгөлөг ялгаварын аргыг ашиглавал

болох ба энд

болон

гишүүд агуулагдсан байна. Энэ тэгшитэглд зөвхөн нэг дууны хурд ашиглагдана. Учир нь D дотор харагдаж байгаа хурднууд нь шинэ хугацааны алхамынх ба эдгээр нь мөн тэгшитгэл 6-аас үүдэн n+1 түвшний даралтаас хамаарах учир энэ тэгшитгэл нь шинэ даралтыг тодорхойлох далд итерици болно. Шингэн агуулж буй нүд бүрийн даралтын өөрчлөлт нь тэгшитгэл 10-ийг хангах ёстой.

Үүнд S нь тэгшитгэл 10-ийн зүүн талыг тэмдэглэсэн ба шинээр итерацилагдаж байгаа даралт нь

гэж тооцоологдох ба нүднүүдийн хажуу тал дээр байрлах хурднууд нь

гэх маягаар шинэчлэгдэнэ. Энд dt, dp гэх мэт нь дифференциал биш бөгөөд өөрчлөлтийг тэмдэглэсэн юм. Омега нь тэгшитгэл 10-ийг тэгшитгэл 13-ын дагуу шийдэхэд бий болж байгаа далд утгыг хянана. Харин чөлөөт гадаргуугын нүднүүдийн хувьд S нь өөр байх ба энд нөгөө шингэний даралт (P_s) болон чөлөөт гадаргуу үүсгэгч шингэний даралт P_n хоорондын градиентыг тооцох хэрэгтэй болно. Энэ схемийг ажиллуулахын тулд функц нь

гэж олдоно. Үүнд, этта=dc/d нь зэргэлдээх нүднүүдийн төв хоорондын зайг интерполяци хийгдэх нүдний төв болон чөлөөт гадаргуугын хоорондын зайтай харьцуулсан харьцаа юм. Ингээд бүрэн итераци нь бүх нүднүүдийг тэгшитгэл 11-ээс 13-аар хурд ба даралтыг засварлаж хийгдэх ба эдгээрт S нь дотоод шингэний нүднүүдэд тэгшитгэл 10-аар, гадаргуугын нүднүүдэд тэгшитгэл 14-өөр өгөгдөнө. Итерацийн нийлэлт нь S-ийн утгын тодорхой бага утга е хүртэл хийгдэнэ. Ерөнхийдөө e=10^-3 байдаг. Зарим тохиолдолд итерацийн нийлэлтийг тэгшитгэл 11-ийн dp –г давуулан тогтворшуулах w утгаар үржүүлж хурдасгадаг. Үүний утга нь 1.8 байх ба 2.0 –оос хэтэрч болохгүй.

Практикт, итерполяцийн утга w нь 1-ээс их болонгуут чөлөөт гадаргуугын нөхцөл тэгшитгэл 14 нь тогтворгүйжилт дагуулдаг учир гадаргуугын нүдэнд интерполяци хийхэд хэрэглэгдэх нүдэн дахь даралтын өөрчлөлтийг доод тогтворшилтоор хангах шаардлагатай байдаг. Ялангуяа нөгөө тогтворшуулагч w нь чөлөөт гадаргуугын хувьд интерполяци хийгдэх нүдэн дээр хэрэглэгдэж байгаа бол

гэж солигдох ёстой болно. Энд этта ба dc нь гадаргуугын нүдэнд хамааралтай ба S нь хөрш нүдний утга юм.

Зураг 3. Чөлөөт гадаргуу дээрх даралтын хязгаарын нөхцлийг тодорхойлоход хэрэглэгдэх хэмжээс, нүдний төрлүүд.

Шингэний эзлэхүүний функцын ойролцоолол

Хугацааны дэвшилд F-ийг тодорхойлох

Шингэний эзлэхүүний функц F нь тэгшитгэл 1-ээр илэрхийлэгдэнэ. Үл шахагдах шингэний хувьд тэгшитгэл 4 нь тэгшитгэл 2-той хослож

гэж илэрхийлэгдэх ба фси= 0 үед r нь 1 байна. Тооцооны нүднүүдийн дагуу тэгшитгэл 16-г нийлбэрчлэх үед нүдэн дэх F-ийн өөрчлөлт нь нүдний ханаар гарах F-ийн цүлхэлтийг тодорхойлно. Өмнө өгүүлсэнчлэн, чөлөөт гадаргуугын тодорхойлолтыг нарийн хадгалахын тулд энэ хэсэгт онцгой анхаарал тавигдах ёстой. Бид донор болон хүлээн авагч гэсэн схемийг хэрэглэнэ. Үндсэн санаа нь доод хэсгийн (доод хашиц) F-ийн тухай мэдээлэл мөн цүлхэлтийн хязгаарын дээд хэсгийн (дээд хашиц) F-ийг харилцах үеийн ойролцоо хэлбэрийг үүсгэхийн тулд хэрэглэх ба дараа нь энэ хэлбэрээ цүлхэлтийн тооцоонд хэрэглэнэ. Шингэний эзлэхүүний техник хэрэглээнд боловсруулж байгаа аргын үндэс нь dt хугацааны алхамын үед нүдний баруун гар талаар F хэмжээтэй шингэн шилжиж байгаа гэдгийг анхаарвал ойлгогдож болох юм. Энэ нүдний талаар гарч байгаа цүлхэлтийн нэгж талбайд ноогдох утга нь Vx=udt ба үүнд u нь нүдний талд нормаль хурд болно. Энэ хурдны тэмдэг нь донор эсвэл хүлээн авагч гэдгийг тогтоох буюу өөрөөр хэлбэл эзлэхүүн алдаж эсвэл авч байгааг илэрхийлнэ. Тухайлбал хэрэв хурд u  нь эерэг бол дээд талын эсвэл зүүн талын нүд нь донор болох ба доод талын эсвэл баруун талын нүд нь хүлээн авагч болно. Нэгж хугацааны алхамд нүдний ханаар шилжиж байгаа цүлхэлтийн хэмжээ нь dF –ийг нүдний хөндлөн талбайгаар үржсэнтэй тэнцүү ба үүнд:

ба

байна. Ганц үсгэн тэмдэглэгээ нь донор болон хүлээн авагчийг, хос тэмдэглэгээ нь аль нь ч байж болохыг заана. Жишээ нь AD нь аль нь ч байж болох ба шингэний урсгалын чиглэлтэй харьцангуй харилцах үеийн байрлалаас хамаарна. Товчоор хэлбэл тэгшитгэл 17-д байгаа минимум нь өгөгдөхөөсөө их хэмжээтэй шингэний цүлхэлт дамжуулагдахаас сэргийлж байхад максимум нь нэмэлт цүлхэлтийг тооцно. Хэрэв хоосрох хэмжээ /1.0-F/ нь байх ёстой хэмжээнээсээ хэтрэх тохиолдолд  максимум чухал үүрэгтэй байна. Донор болон хүлээн авагч нүднүүд нь зураг 5а-д дүрслэгдсэн байгаа. Ad=D үед цүлхэлт нь ердийн донор нүдний утгатай

тэнцүү байх ба энд донор нүдний F-ийн утга нь F –ийн нэвтэрч буй нүдний хажуу хананы фракцийн талбайг тодорхойлоход хэрэглэгдэнэ /зураг 5б/. Тооцооны тогворшилт нь /Vx/ нь dx –ээс бага байхыг шаардах ба ингэснээр энэ тохиолдолд донор нүд хоосрох боломжгүй юм.

Зураг 5. Шингэний эзлэхүүн F –ийн шилжилтэнд хэрэглэгдэх чөлөөт гадаргуугын хэлбэрүүд. Донор ба хүлээн авагч схем зургийн а-д харагдах ба энд тасархай зураас нь нийт эзлэхүүний зүүн хязгаар нь шилжсэн байгааг илтгэнэ. зургийн b-d үзүүлсэн хэсэгт хэрээслэсэн хэсэг нь шилжих гэж буй бодит F –ийн хэмжээг харуулна.

AD=A байх үед хүлээн авагчид байх F-ийн утга нь F-ийн урсан орж ирж байгаа хажуу хананы фракцын талбайг тодорхойлоход хэрэглэгдэнэ. Энэ тохиолдолд (зураг 5с) донор нүдэнд байгаа бүх F шингэн шилжинэ. Учир нь тасархай зураас болон цүлхэж буй хана хооронд байгаа бүх зүйл хүлээн авагч нүдрүү нүүнэ. Энэ бол тэгшитгэл 17 байгаа минимумын жишээ юм. Зураг 5-ийн d тохиолдолд Fa/Vx/-аас илүү их хэмжээний F шингэн шилжсэн байх ёстой. Энэ бол максимумын жишээ юм. Тухайлбал тасархай зураас болон цүлхэлтийн хязгаар хооронд байх илүүдэл шингэн нь тэгшитгэл 17 дахь CF утгатай тэнцүү байна. Хүлээн авагч болон донор нүднүүдийн аль нь F –ийн урсгал дахь фракцын тайлбайг тодорхойлогдоход хэрэглэгдэх эсэх нь гол гадаргуугын чиглэлээс хамаарна. Хэрэв гударгуу нь ихэвчлэн өөртэйгээ нормалаар шилжиж байвал хүлээн авагч нүд хэрэглэгдэнэ. Эсрэг тохиолдолд донор нүд  хэрэглэгдэнэ. Гэхдээ хүлээн авагч нүд нь хоосон бол эсвэл донор нүдний дээд хэсгийн нүд хоосон бол хүлээн авагч нүдний F утга нь гадаргуугын чиглэлээс үл хамааран цүлхэлтийн утгыг тодорхойлох болно. Энэ нь донор нүд ямар нэгэн шингэний эзлэхүүн доод хашиц дахь нүдрүү орохоос өмнө дүүрсэн байх ёстой гэсэн үг юм. Гадаргуугын чиглэл дээр хийгдэх тестийн шалтгаан нь хэрэв хүлээн авагч нүд нь цүлхэлтийг тодорхойлоход үргэлж хэрэглэгддэг бол гадаргуугын буруу налуу долгион үүсдэгтэй холбоотой юм. Жишээ нь жижиг долгионтой хэвтээ гадаргуу эерэг х чиглэлд хөдөлж байгаа гэж үзвэл доод талд /хүлээн авагч/ байрлах цүлхэлт дээр суурилсан F –ийн утга нь алхамын тасралтруу долгионыг огцом болгож өгнө. Энэ нөлөөлөлд хүлээн авагч нүд нь тогтворгүйжилт бий болгох ба энэ нь F –ийн сөрөг утгатай диффузыг бий болгоно. Өөрөөр хэлбэл сөрөг коэффициенттэй диффузи маягийн шилжилт бий болно. Энэ тогтворгүйжилт хязгааргүйгээр өсөхгүй, гэсэнч мин болон макс тест цүлхэлтийн тодорхойлолтонд хэрэглэгдэнэ. Нэгэнт цүлхэлт дээрхээр тодорхойлогдчихвол хязгаарын талбайгаар үржигдэж донор нүднээс F хэмжээтэй шингэн хасагдаж, хүлээн авагч нүдрүү F хэмжээтэй шингэн нэмэгдэж орно.

Балансын тохируулга

Дээрх аргаар тодорхойлогдсон F –ийн шинэ утга нь ямар нэгэн байдлаар 0-ээс бага эсвэл 1-ээс их байх тохиолдол байна. Тийм учраас шилжилтийн тооцоолол хийгдсэний дараагаар нэг удаагийн явалт бүх нүдэн дээгүүр явагдаж 0-ээс бага утгыг 0 болгож, 1-ээс их утгыг нэг болгоно. Тооцооллын үеэр хийгдэх эдгээр тохируулгаар гарч ирсэн шингэний эзлэхүүний хуримтлагдсан өөрчлөлт нь хугацааны алхам бүрт хадгалагдаж, хэвлэгдэж байх хэрэгтэй. Шингэний эзлэхүүн F –д өөр нэг тохируулга хэрэгтэй, ингэснээр хязгаарын нүдний тугаар хэрэглэгдэж болох юм. Өөрөөр хэлбэл хэрэв F –ийн утга нь эфсо_F утгаас бага болвол хоосон нүд болох ба F –ийн утга нь F-эфсо_F утгаас их болвол шингэнээр дүүрсэн нүд болно. Хэрэв шилжилтийн дараагаар нүд нь эфсо_F –ээс бага утгатай F –тэй бол энэ F-ийг 0 болгож бусад хөрш шингэний нүднүүдийн F утгыг 1.1эфсо_F –ээр багасгаж гадаргуугын нүднүүд болгоно. Эзлэхүүний утга F дахь эдгээр өөрчлөлт нь хуримтлагдсан эзлэхүүний өөрчлөлтөнд багтсан байх хэрэгтэй. Зуун циклийн дараах эзлэхүүний алдаа нь нийт F эзлэхүүний зууны 1 хувьиас хэтрэхгүй бол сайн.  

Донор-хүлээн авагчийн схемийг Фортран хэлээр кодолсон байна.

© Hosoyamada Tokuzo

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

do while(t.lt.supt) ! time advancing ilop=ilop+1 do i=1,isize ! f1 for temporal change f0 for previous do j=1,isize f1(i,j)=f0(i,j) enddo enddo do i=2,isize-1 do j=2,isize-1 yi(i,j)=(f0(i ,j-1)+f0(i,j)+f0(i ,j+1))*h !dundaj*h xi(i,j)=(f0(i-1,j )+f0(i,j)+f0(i+1,j ))*h enddo enddo do i=3,isize-2 do j=3,isize-2 sl1=abs(2.0*(yi(i+1,j )-yi(i-1,j ))/(4.0*h)) sl2=abs(2.0*(xi(i ,j+1)-xi(i ,j-1))/(4.0*h)) if(sl1.le.sl2) then ihv(i,j)=0 !horizontal else ihv(i,j)=1 !vertical endif enddo enddo do i=3,isize-3 do j=3,isize-3 ! X-direction dv=abs(u(i,j)*dt) if(u(i,j).ge.0) then id=i ! donner ia=i+1 ! acceptor else id=i+1 ! donner ia=i ! acceptor endif if(ihv(id,j).eq.1) then !donner-n huvid bosoo iad=ia !a-d n accep bolno else !hewtee iad=id !a-d n donner bolno endif if(f0(ia,j).lt.emikro) iad=id !a-d n donner bolno cf=dmax1( (1.0-f0(iad,j))*dv-(1.0-f0(id,j))*h,0. ) ! text df=dmin1( f0(iad,j) *dv+cf ,f0(id,j)*h ) ! text f1(id,j)=f1(id,j)-df/h !donner-s hasna f1(ia,j)=f1(ia,j)+df/h !accep-t nemne ! y-direction dv=abs(v(i,j)*dt) !flux utga if(v(i,j).ge.0) then jd=j !j donner ja=j+1 !j acceptor else jd=j+1 ja=j endif if(ihv(i,jd).eq.1) then !bosoo jad=jd else jad=ja endif if(f0(i,ja).lt.emikro) jad=jd cf=dmax1((1.0-f0(i,jad))*dv-(1.0-f0(i,jd))*h, 0. ) df=dmin1( f0(i,jad) *dv+cf ,f0(i,jd)*h ) f1(i,jd)=f1(i,jd)-df/h f1(i,ja)=f1(i,ja)+df/h enddo enddo do i=1,isize do j=1,isize if(f1(i,j).lt.emikro) then f1(i,j)=0 if(f1(i-1,j ).ge.1.0) f1(i-1,j )=f1(i-1,j )-1.1*emikro if(f1(i+1,j ).ge.1.0) f1(i+1,j )=f1(i+1,j )-1.1*emikro if(f1(i ,j-1).ge.1.0) f1(i ,j-1)=f1(i ,j-1)-1.1*emikro if(f1(i ,j+1).ge.1.0) f1(i ,j+1)=f1(i ,j+1)-1.1*emikro endif f0(i,j)=dmin1(1.0,dmax1(f1(i,j),0.)) enddo enddo t= t+dt it=it+1 enddo ! temporal time advancing

Шингэний эзлэхүүний шилжилтийг илэрхийлэх адвекцийн тэгшитгэлийг шийдэх аргагчлалыг хайж байна.