Шингэний чөлөөт гадаргууг ЛБА тооцоолох
Усны барилгын хувьд шингэн ихэнхдээ агаартай харьцах ба шингэн ба агаарын харилцах үеийг чөлөөт гадаргуу гэж нэрлэнэ. Хэлхээст тооцооны аргад уг чөлөөт гадаргууг тодорхойлоход төвөгтэй байх ба чөлөөт аргад ямар нэгэн нэмэлт техник ашиглалгүйгээр ойролцоолох боломжтой байдаг. ЛБА нь мөн л хэлхээст аргын ангилалд багтах учир нэмэлт загвар ашиглаж усны чөлөөт гадаргууг нэг фазат байдлаар тооцох боломжтой. Төгсгөлөг аргуудын хувьд шингэний эзлэхүүний арга эсвэл шийдлийн алгоритмоор шингэний чөлөөт гадаргууг ойролцоолох боломжтой.
Шингэний цалгилт
Усан баганын задралын шинжилгээ
Шингэний чөлөөт гадаргууг тооцоолох
Чөлөөт гадаргуутай шингэний симуляцийг нэг
фазат ЛБА-аар тооцоход бодлогийн хүрээг шингэн болон хий агуулсан ялгаатай
бүсүүдтэй байхыг шаарддаг. Энэ нь кодонд тухайн нүд нь зөвхөн хий агуулах бол
хийн нүд гэсэн тэмдэглэгээг (туг) хийхийг шаардана. Хатуу
гадаргуугын нүдтэй ижил байдлаар эдгээр хийн нүднүүдийн түгэлтийн функц нь
симуляцийн үед тооцоонд оролгүй орхигдоно. Гэхдээ хатуу гадаргуугын нүднээс ялгаатай
тал нь симуляцид зарим нүдэн дахь шингэн нь уг хоосон бүсрүү хөдлөх боломжтой. Шингэний
энэ хөдөлгөөнийг тодорхойлохын тулд өөр нэг нүдийг авч үзэх шаардлагатай. Энэ
бол харилцах үеийн нүд юм. Эдгээр нүд нь хаагдсан нэгж зузаантай давхрага байх
ба ямагт шингэн ба хийн нүдээр хоёр талдаа хүрээлэгдэнэ. Шингэний чөлөөт
гадаргуугийн хөдөлгөөнийг тооцоолохын тулд ЛБА-д нэмэлтээр дараах гурван
алхамын тооцоог хийх хэрэгтэй. Үүнд: харилцах үеийн хөдөлгөөн, шингэний
харилцах үе дээрх хязгаарын нөхцөл, тооцооны нүдний төрлийг дахин тогтоох гэх
мэт болно. Дараах зураг нь тооцооны гол алгоритмыг дүрслэн харуулсан байна.
4.1 Чөлөөт гадаргуугын хөдөлгөөн
Шингэний чөлөөт гадаргуугын хөдөлгөөн нь
массын шилжилтээр тодорхойлогдох ба уг масс нь тооцооны нүд бүр дээр тодорхой
утга авна. Массын шилжилтийг тооцоолоход хоёр нэмэлт хэмжигдэхүүнийг тооцооны
нүд бүр дээр тодорхойлох шаардлага гарах ба эдгээрийг массын хэмжээ m болон
шингэний фракц e гэж нэрлэнэ. Шингэний эзлэхүүний фракц нь тухайн нүдний масс
ба шингэний нягтын харьцаагаар тодорхойлогдоно.
Шингэний эзлэхүүний арга шиг харилцах
үеийн хөдөлгөөн нь нүд хоорондын массын цүлхэлтийн хэмжээгээр тооцоологдоно.
Гэхдээ тодорхой тооны эгэл хэсгийг багтаасан түгэлтийн функц байгаа учир массын
өөрчлөлт нь энэхүү загварт хөрш хоёр нүдний хооронд урсах түгэлтийн функцуудын
ялгавараар тодорхойлогдоно. Хэрэв x+deltat*ci зүгт шингэний нүд байх бөгөөд идэвхитэй байгаа нүд нь харилцах үеийн
нүд бол массын өөрчлөлт нь дараах байдлаар тодорхойлогдоно.
Эхний түгэлтийн функц нь энэ хугацааны
алхамд идэвхит байгаа нүдрүү нүүж ирэх шингэний массын хэмжээг илтгэх ба хоёр
дахь функц нь энэ нүдийг орхиж байгаа шингэний массын хэмжээг илэрхийлнэ. Харин
хоёр харилцах нүд хоорондын масс солилцоо нь тус бүрийн шингэнээр дүүргэгдсэн
талбайн хамаарлаар шилжих болно. Шингэнээр дүүрсэн хэсгийн талбайг хоёр нүдний
хувьд шингэний фракцуудын утгыг дундажилж тодорхойлж болно. Тиймээс тэгшитгэл
4.2 нь дараах байдалтай болно.
Дээрх тэгшитгэлд байх массын цүлхэлтийг
харилцах үеийн залруулга хэсэгт өгөгдсөн тодорхой дүрмийн дагуу (хүснэгт) тодорхойлно. Тэгшитгэл 4.2 болон 4.3 – ийг харвал
нүдрүү ирж байгаа болон нүдийг орхиж байгаа массын харьцаа нь тэгш хэмтэй болох
нь харагдаж байна. Иймд масс хадгалагдах хууль дараах байдлаар батлагдана.
Хөрш шингэний нүдтэй харилцах үеийн нүдний
хувьд массын солилцоо нь урсгалын үе дахь түгэлтийн функцын солилцоотой тохирох
хэрэгтэй ба шингэний нүдний хувьд нэмэлт тооцоолол шаардлагагүй. Тэдгээрийн
шингэний фракцууд нь үргэлж нэгтэй тэнцүү байх ба масс нь өөрсдийнх нь нягттай
тэнцүү байх учир массын өөрчлөлт гарахгүй. Бүх чиглэл дахь массын өөрчлөлтийн
утга нь харилцах үеийн нүдний масс дээр нэмэгдэх ба гарсан массын утга нь
дараачийн хугацааны алхамын тооцоонд хэрэглэгдэнэ.
4.2 Чөлөөт гадаргуугын хязгаарын нөхцөл
Дээр тайлбарласанчлан хоосон нүдэнд байх
түгэлтийн функцууд хэзээ ч тооцоонд орохгүй. Гэсэн ч харилцах үеийн нүдэнд ядаж
нэг хоосон нүдтэй хөрш байх учир тухайн нүднээс нүүж ирэх түгэлтийн функцыг
тооцох шаардлагатай. Тиймээс урсах алхамд шингэн ба харилцах үеээс түгэлтийн
функцууд энгийн байдлаар нүүж ирнэ. Харин хоосон буюу хийн нүднээс ирж байгаа
түгэлтийн функцууд нь хязгаарын нөхцлөөр дахин бүрэлдэх ёстой. Энэ хязгаарын
нөхцөл нь нэмэлт нүдний үе гэх мэтийг шаардахгүй идэвхит зангилаан /нүд/ дээрээ
тооцогдоно. Атмосферийн даралтыг тооцохын тулд агаарын нягтыг нормаль байдлаар rhoa=1 гэж өгөх ба
эдгээр нь мөн лавлах нягт болон даралт болох болно. Цаашилбал зунгааралт нь
хийн хувьд маш их байхад шингэний хувьд маш бага байна. Энэ нь шингэний фазад
тооцооллийг хийж хийн фазыг алгасах нөхцөлийг бүрдүүлж байгаа юм. Түгэлтийн
функцын хувьд хэрэв x+ci нь хийн нүд байвал
гэж
тодорхойлогдоно. Үүнд u нь t хугацааны x байрлал дахь урсгалын хурд юм. Энэ нь
түгэлтийн функцын 1 - р эрэмбийн моментоор илэрхийлэгдэж тэндээс олдоно.
Шингэний харилцах
үед үйлчилж байгаа атмосферийн даралт нь агаарын нягтыг хэрэглэсэн тэнцвэрт
түгэлтийн функцаар тодорхойлогдоно. Хоосон нүдтэй бүх хөршийн чиглэлд тэгшитгэл
4.5 -ыг хэрэглэж харилцах үеийн нүдэнд түгэлтийн функцын бүрдлийг бий болгоно. Гэхдээ
харилцах үеийн тал бүр дээр байх балансыг хадгалахын тулд харилцах үеийн нормал
чиглэлээс ирж байгаа түгэлтийн функцууд мөн дахин байгуулагдах ёстой. Тиймээс
хэрэв хоосон болон дараах нөхцөл биелсэн нүднээс нүүж ирэх түгэлтийн функцуудыг
хязгаарын нөхцлийн тусламжтайгаар байгуулна.
Энд j,k,l нь торны
хэвтээ ба босоо тэнхлэгүүдийн дугаарыг заана. Эндээс нормаль нь чиглэл бүр дахь
шингэний фракцийн төвийн ялгавараар ойролцоологдож байна.
Одоо харилцах үе
дээрх бүх түгэлтийн функцууд байгуулагдсан учир мөргөлдөөнийг явуулах бүрэн
боломжтой. Мөргөлдөөний үед хурд ба нягт тооцогдох учир дараах харьцааг
шалгахад боломжтой болж байна.
Энд нэмэлт зайн шилжүүлэг
ашиглах нь шинээр бий болсон харилцах үеийн нүднүүдийг дараачийн ЛБ-ны алхамд
дахин хувирхаас сэргийлэх ба 0 болон 1 гэсэн утгын орон k=10-3 гэсэн босгыг
хэрэглэж байгаа юм. Хоосон болон дүүрсэн нүднүүдийн гэнэтийн хувирлын оронд
тэдгээрийн байрлал нь ямар нэгэн түр жагсаалтанд хадгалагдаж (нэг нь хоосон
нүдэнд, нөгөөх нь дүүрсэн нүдэнд) бүхэл процесс бүх нүдийг хамарсаны дараа жинхэнэ
шилжүүлэг хийгдвэл зохино.
4.3 Нүднүүдийн төрлийг тогтоох, тооцоонд дахин
бэлтгэх
Энэ алхам нь бүх
нүднүүд шинэчлэгдсэн үед дараах хоёр шинж чанарыг батлахын тулд хийгдэнэ. Үүнд
харилцах үеийн давхрага нь хаагдсан байх ёстой. Нэгэнт дүүрсэн болон хоосорсон
харилцах үеийн нүд нь тэдгээрийн хувирах ёстой нүдний төрөврүү хувирах
хэрэгтэй. Дээрээс нь нэмээд массын илүүдэл нь хувирсан нүднүүдээс шилжүүлгийн
үед бусад харилцах үеийн нүднүүдрүү тараагдах ёстой. Тэгшитгэл 4.7-оор массыг
жишихэд нүд дүүрсэн ч, хоосорсон ч ямар нэг байдлаар илүүдэл масстай байна. Уг
масс нь эерэг ба сөрөг гарах ба хөрш харилцах үеийн нүднүүдрүүгээ тараагдах
хэрэгтэй юм.
Эхлээд дүүрсэн
нүдний бүх хөрш нүднүүд бэлтгэгдсэн байна. Бүх хөрш хоосон нүднүүд харилцах
үеийн нүд болж хувирна. Идэвхит нүдийг хүрээлж буй харилцах үеийн болон
шингэний нүднүүдээс дундаж нягт rhoav болон дундаж
хурднууд uav тодорхойлогдсон
байна. Хоосон нүдний түгэлтийн функцууд нь тэнцвэрт түгэлтийн функын
томъёоллоор fieq(rhoav,uav) гэж шинээр
байгуулагдана. Энд зарим нэг хоосорсон гэсэн жагсаалтанд байгаа харилцах үеийн
нүд нь хувирч буй шингэний нүдэнд хөрш болох учраас харилцах үеийн нүд хэвээр
үлдэж түр жагсаалтаас гарах хэрэгтэй болно. Ижил хугацааны үед дүүрсэн нүдний түр
туг нь шингэний туг болж өөрчлөгдөнө. Үүнтэй адилаар хоосорсон нүдний бүх хөрш
шингэний нүд нь харилцах үеийн нүдрүү шилжиж, байгаа түгэлтийн функцууд нь уг
шинэ харилцах үеийн нүдэнд хэрэглэглэх болно. Ингээд хоосорсон харилцах үеийн
нүд нь одоо жинхэнэ хоосон нүд гэж тэмдэглэгдэх боломжтой боллоо. Хоёр дахь
ажил нь илүүдэл масс mex нь хүрээлж буй харилцах үеийн нүднүүдрүү
хуваарилагдах ёстой ажил юм. Илүүдэл масс нь хоосорч буй нүдний хувьд түүний
масстай, харин дүүрч буй нүдний хувьд түүний масс ба нягтын ялгавартай тэнцүү
байна.
Хоосорсон харилцах
үеийн нүдний масс сөрөг байна гэдэг нь дүүрсэн нэгний хувьд масс нь нягтаасаа
их байгаатай адил шингэний харилцах үе одоогийн нүднээс цааш ухарч шилжиж байна
гэсэн үг юм. Үүнийг тооцохын тулд илүүдэл масс нь хөрш харилцах үерүүгээ тэгш
хуваарилагдахгүй гэхдээ харилцах үеийн нормалын чиглэлийн дагуу жинлэгдэх ёстой
болно.
Энд эта нийт нь
чиглэл бүр дахь эта-гийн нийлбэртэй тэнцүү бөгөөд эта нь дараах байдлаар
тодорхойлно.
Хажуугийн харилцах
үеийн массын өөрчлөлт шиг нүдний шингэний
фракц нь мөн л өөрчлөгдөх хэрэгтэй. Иймд масс тараалтын дараагаар шингэний
эзлэхүүний фракцыг тодорхойлох хэрэгтэй. Энэ хүртэл тайлбарлагдсан алхамуудын
хувьд дүүрсэн ба хоосорсон нүдний хувиралын дарааллаас үл хамааран ижил үр дүнг
гаргаж авах боломжтой. Тэгэхээр хоосон нүднүүдийн интерполяци нь өөрсдөө шинэ бус
харилцах үеийн нүднүүдийн утгын хувьд хийгдэж байж болох юм. Нэгэнт нүдний
шилжүүлгүүд хийгдсэн бол идэвхит тор нь хүчин төгөлдөр болж бүх нүдний хувьд
шинэ гогцоо давтамж эхлэх болно.
Харилцах үеийн нүдний залруулга
Дээр дурьдагдсан
алгоримт нь чөлөөт гадаргуугын хөдөлгөөнийг илэрхийлэхэд хүрэлцээтэй юм. Гэхдээ
шингэн хөдөлж байхад ганц харилцах үеийн нүд үлдэх эсвэл харилцах үеийн нүд
шингэн дотор хүрээлэгдэх гэх мэт нь загварын согогыг бий болгоно. Эдгээр нь
шингэний симуляцийг тасалдуулахгүй боловч согог мэт харагдсаар байна. Үүнийг
намжаахын тулд дараах дүрэм алгоритмд хэрэглэгдэх болно. Үндсэн санаа нь хөрш
шингэний нүдийг хоосон болгохгүйгээр харилцах үеийн нүдийг хэвээр үлдээх эсвэл
хоосон хөршүүдийг дүүргэхгүйгээр харилцах үеийн нүдийг хэвээр байлгах гэх мэт
юм. Энэ нь f-i(x+deltat*ci,t)-fi(x,t) гэсэн нэмэлтийг
тэгшитгэл 5-д хүснэгт 1-ийн дагуу хийснээр зохицуулагдана. Ховор тохиолдолд
харилцах үеийн нүднүүд нь хэдий харилцах үе байсан ч энгийн байдлаар дүүрсэн
эсвэл хоосорсон нүд болж хувирах ёстой. Тиймээс хэрэв нүд нь ямарч шингэн
хөршгүй бөгөөд үүний масс нь 0.1ро – оос доош орсон бол энэ хугацааны алхамд
хоосорсон нүд гэж тооцогдоно. Харин нүд нь ямарч шингэн хөршгүй бөгөөд масс нь
0.9ро-оос дээш байвал дүүрсэн гэж тооцогдоно.
Дээрх алгоритм болон тэгшитгэл томъёоллыг хэрхэн кодлох талаар авч үзье. Ерөнхий ЛБа-ын код дээр массын өөрчлөлт, нүдний төрөл түүний өөрчлөлт, хязгаарын нөхцөл гэх мэт дэд програмууд нэмэгдэнэ.
Массын өөрчлөлтийг тодорхойлох дэд програм
Массын өөрчлөлтийг тодорхойлох дэд програм
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 | do i=0,n !range on domain or range inside the domain do j=0,m asum=0.0 if(flag(i,j).eq.1) then ! fluid node !1 ma(i,j)=rho(i,j) do k=1,8 op=opp(k) if(flag(i+cx(k),j+cy(k)).eq.2.or. & flag(i+cx(k),j+cy(k)).eq.1) then delm(k,i,j)=f(op,i+cx(k),j+cy(k))-f(k,i,j) endif asum=asum+delm(k,i,j) enddo ma(i,j)=ma(i,j)+asum else if(flag(i,j).eq.2) then ! interface node do k=1,8 op=opp(k) if(flag(i+cx(k),j+cy(k)).eq.3) then delm(k,i,j)=0. else if(flag(i+cx(k),j+cy(k)).eq.1) then delm(k,i,j)=f(op,i+cx(k),j+cy(k))-f(k,i,j) else if(flag(i+cx(k),j+cy(k)).eq.2) then ! mass exchange between cells nearest IF call unwant(i,j,k,delm) ! think delm(k,i,j)=0.5*(e(i+cx(k),j+cy(k))+e(i,j)) & *delm(k,i,j) ! print*, 'unwanted subroutine is used',i,j,k,delm(k,i,j) endif asum=asum+delm(k,i,j) enddo ! print*, 'asum', asum ! endif ma(i,j)=ma(i,j)+asum if(isnan(ma(i,j))) stop ! write(*,*) 'mass and density', ma(50,40), rho(50,40) endif enddo enddo |
Дээрх кодны хэсэгт flag нь нүдний төрлийг заах ба 1-ус, 2-харилцах үе, 3-агаар байна.
Энд байх unwant /25-р мөр/ хэмээх дэд програмын үндсэн үйлдлийг доорхоор харуулж байна.
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 | op=opp(k) if(flagneig(i,j).eq.2) then ! second line ! do k=1,8 if(flagneig(i+cx(k),j+cy(k)).eq.2) then delm(k,i,j)=f(op,i+cx(k),j+cy(k))-f(k,i,j) else ! flagneig is 1 and 3 delm(k,i,j)=-f(k,i,j) endif print*, 'I am second line' ! enddo else if(flagneig(i,j).eq.3) then ! third line if(flagneig(i+cx(k),j+cy(k)).eq.3) then delm(k,i,j)=f(op,i+cx(k),j+cy(k))-f(k,i,j) else ! flag is 1 and 2 delm(k,i,j)=f(op,i+cx(k),j+cy(k)) endif print*, 'I am third line' else if(flagneig(i,j).eq.1) then ! first line if(flagneig(i+cx(k),j+cy(k)).eq.2) then delm(k,i,j)=f(op,i+cx(k),j+cy(k)) else if(flagneig(i+cx(k),j+cy(k)).eq.3) then delm(k,i,j)=-f(k,i,j) else if(flagneig(i+cx(k),j+cy(k)).ne.0 & .or.flagneig(i+cx(k),j+cy(k)).eq.1) then delm(k,i,j)=f(op,i+cx(k),j+cy(k))-f(k,i,j) endif endif |
Энэхүү дэд програмд байх flagneig нь тухайн нүдний хөршийн талаарх мэдээллийг агуулах ба ямарч шингэн нүдний хөрш байхгүй бол 2, ямарч хоосон нүдний хөрш байхгүй бол 3, шингэн болон хоосон хөрштэй мөн харилцах үеийн хөрштэй харилцах үеийн нүд бол стандарт 1 гэсэн утгыг авна.
Массын өөрчлөлтийг тооцсоны дараа ямар харилцах нүд шингэн эсвэл агаарын нүд болохыг шийдэж илүүдэл массыг тараах, шингэний фракцыг тооцох, нүдний төрлийг дахин тогтоох гэх мэт чухал үйлдлүүдийг агуулсан дэд програмыг доорхоор харуулав.
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 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 | ! temporary flag setting do j=0,m do i=0,n if(flag(i,j).eq.2) then if(ma(i,j).gt.1.003*rho(i,j)) then tempflag(i,j)=4 ! IF to F /IFF/ print*, 'IFFa cell i=', i, 'j=', j else if(ma(i,j).lt.-0.003*rho(i,j)) then tempflag(i,j)=5 ! IF to G /IFG/ print*, 'IFGb cell i=', i, 'j=', j endif if(ma(i,j).ge.(1.-0.1)*rho(i,j).and.flagneig(i,j).eq.2) then tempflag(i,j)=4 !IF to F print*, 'IFFc cell (lonelytresh) i=', i, 'j=', j else if(ma(i,j).lt.0.1*rho(i,j).and.flagneig(i,j).eq.2) then tempflag(i,j)=5 !IF to G print*, 'IFGd cell (lonelytresh 1) i=', i, 'j=', j else if(flag(i,j).eq.2.and.flagneig(i,j).eq.2) then tempflag(i,j)=5 !IF to G AND RECONSIDER IT!!!!! print*, 'IFGe cell (lonelytresh 2) i=', i, 'j=', j endif endif enddo enddo ! changing neighboring flags bacause of ordering current cell flag change do j=0,m ! rethink do i=0,n if(tempflag(i,j).eq.4) then ! IFF do k=1,8 if(flag(i+cx(k),j+cy(k)).eq.3) then ! calculate surrounding rho, u and then f f(k,i+cx(k),j+cy(k))=w(k)*rho(i+cx(k),j+cy(k)) call reinit(i+cx(k),j+cy(k),k) !problem is here ! f(k,i+cx(k),j+cy(k))=f(k,i+cx(k),j+cy(k)) flag(i+cx(k),j+cy(k))=2 print*, 'new cells are reinitiialized' endif enddo do k=1,8 if(tempflag(i+cx(k),j+cy(k)).eq.5) then flag(i+cx(k),j+cy(k))=2 tempflag(i+cx(k),j+cy(k))=0 ! no longer tempflag print*, 'neighborhood for filled cell unlisted from IFG' endif enddo endif enddo enddo do j=0,m do i=0,n if(tempflag(i,j).eq.5) then ! IFG do k=1,8 if(flag(i+cx(k),j+cy(k)).eq.1) then flag(i+cx(k),j+cy(k))=2 print*, 'IF cell changed to Gas cell' endif enddo do k=1,8 if(tempflag(i+cx(k),j+cy(k)).eq.4) then flag(i+cx(k),j+cy(k))=2 ! last addition tempflag(i+cx(k),j+cy(k))=0 ! no longer tempflag print*, 'neighborhood for emptied cell unlisted from IFF' endif enddo endif enddo enddo ! calculation of sum of weight. do i=0,n ! range remind do j=0,m ifnum(i,j)=0 ! if(flag(i,j).eq.2) then ! interface node sumw(i,j)=0.0 ! sum of weight if(tempflag(i,j).eq.4) then do k=1,8 if(flag(i+cx(k),j+cy(k)).eq.2) ifnum(i,j)=ifnum(i,j)+1 if(nx(i,j)*cx(k)+ny(i,j)*cy(k).gt.0.) then weight=nx(i,j)*cx(k)+ny(i,j)*cy(k) else weight=0. endif sumw(i,j)=sumw(i,j)+weight enddo ! sumw(i,j)=sumv print*, 'total weight of IFF', sumw(i,j) endif if(tempflag(i,j).eq.5) then do k=1,8 if(flag(i+cx(k),j+cy(k)).eq.2) ifnum(i,j)=ifnum(i,j)+1 if(nx(i,j)*cx(k)+ny(i,j)*cy(k).lt.0.) then weight=-(nx(i,j)*cx(k)+ny(i,j)*cy(k)) else weight=0. endif sumw(i,j)=sumw(i,j)+weight enddo ! sumw(i,j)=sumv print*, 'total weight of IFG', sumw(i,j) endif ! endif enddo enddo ! excess mass determination and distribution do j=0,m ! range remind do i=0,n if(tempflag(i,j).eq.4) then !IFF mex(i,j)=ma(i,j)-rho(i,j) ma(i,j)=rho(i,j) do k=1,8 if(sumw(i,j).gt.0.) then if(nx(i,j)*cx(k)+ny(i,j)*cy(k).gt.0.) then !if 3 wei=nx(i,j)*cx(k)+ny(i,j)*cy(k) else !if 3 wei=0.0 endif if(flag(i+cx(k),j+cy(k)).eq.2) ma(i+cx(k),j+cy(k)) & =ma(i+cx(k),j+cy(k)) & +mex(i,j)*(wei/sumw(i,j)) !rethink after development ! fm(k,i,j)=mex(i,j)*(wei/sumw(i,j)) print*, 'eta_total>0 mass distributed IFF' else if(ifnum(i,j).gt.0) then if(flag(i+cx(k),j+cy(k)).eq.2) ma(i+cx(k),j+cy(k)) & =ma(i+cx(k),j+cy(k))+mex(i,j)/float(ifnum(i,j)) ! fm(k,i,j)=mex(i,j)*(wei/sumw(i,j)) print*, 'ifnum>0 mass distributed IFF' endif enddo else if(tempflag(i,j).eq.5) then ! IFG mex(i,j)=ma(i,j) ma(i,j)=0. do k=1,8 if(sumw(i,j).gt.0.) then if(nx(i,j)*cx(k)+ny(i,j)*cy(k).lt.0.) then wei=-(nx(i,j)*cx(k)+ny(i,j)*cy(k)) else wei=0.0 endif if(flag(i+cx(k),j+cy(k)).eq.2) ma(i+cx(k),j+cy(k)) & =ma(i+cx(k),j+cy(k)) & +mex(i,j)*(wei/sumw(i,j)) ! fm(k,i,j)=mex(i,j)*(wei/sumw(i,j)) print*, 'eta_total>0 mass distributed IFF' else if(ifnum(i,j).gt.0) then if(flag(i+cx(k),j+cy(k)).eq.2) ma(i+cx(k),j+cy(k)) & =ma(i+cx(k),j+cy(k))+mex(i,j)/float(ifnum(i,j)) ! fm(k,i,j)=mex(i,j)*(wei/sumw(i,j)) print*, 'ifnum>0 mass distributed IFF' endif enddo endif ! interfaces enddo enddo ! mass determination on cells do j=0,m do i=0,n ! if(flag(i,j).eq.2) then do k=1,8 ! ma(i+cx(k),j+cy(k))=ma(i+cx(k),j+cy(k))+fm(k,i,j) if(tempflag(i,j).eq.4) then flag(i,j)=1 tempflag(i,j)=0 ! if(flag(i-cx(k),j-cy(k)).eq.3) flag(i-cx(k),j-cy(k))=2 print*, 'temporary flag changed to real flag IFF' else if(tempflag(i,j).eq.5) then flag(i,j)=3 tempflag(i,j)=0 ! if(flag(i-cx(k),j-cy(k)).eq.1) flag(i-cx(k),j-cy(k))=2 print*, 'temporary flag changed to real flag IFG' endif enddo enddo enddo ! flag change do j=0,m do i=0,n if(flag(i,j).eq.2) then jF=0 do k=1,8 if(flag(i-cx(k),j-cy(k)).ne.1) jF=jF+1 enddo if(jF.eq.8) flag(i,j)=3 endif enddo enddo ! volume fraction change do j=0,m do i=0,m if(flag(i,j).eq.1) then e(i,j)=1. else if(flag(i,j).eq.2) then if(rho(i,j).gt.0) then e(i,j)=ma(i,j)/rho(i,j) if(e(i,j).ge.1.) then e(i,j)=1.0 else if(e(i,j).lt.0) then e(i,j)=0.0 endif else e(i,j)=0.5 endif else if(flag(i,j).eq.3) then e(i,j)=0. endif enddo enddo ! set neigboring information on IF cell do j=1,m-1 do i=1,n-1 iamF=0 iamIF=0 iamG=0 do k=1,8 if(flag(i-cx(k),j-cy(k)).eq.1) then iamF=iamF+1 else if(flag(i-cx(k),j-cy(k)).eq.3) then iamG=iamG+1 else if(flag(i-cx(k),j-cy(k)).eq.2) then iamIF=iamIF+1 endif enddo if(iamF.ge.1.and.iamG.ge.1.and.iamIF.ge.1) then flagneig(i,j)=1 else if(iamF.ge.1) then flagneig(i,j)=3 else if(iamG.ge.1) then flagneig(i,j)=2 ! else if(flag(i,j).eq.2.and.iamIF.lt.2) then ! flagneig(i,j)=4 endif enddo enddo print*, 'cells state are changed' |
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 | ! sum of density and velocity anum=0. rhave=0. averh=0. uuave=0. aveuu=0. uvave=0. aveuv=0. ii=i-cx(ki) jj=j-cy(ki) do k=1,8 if(flag(ii+cx(k),jj+cy(k)).ne.3) then anum=anum+1. rhave=rhave+rho(ii+cx(k),jj+cy(k)) uuave=uuave+u(ii+cx(k),jj+cy(k)) uvave=uvave+v(ii+cx(k),jj+cy(k)) endif enddo averh=rhave/anum ! average density of neighbors aveuu=uuave/anum aveuv=uvave/anum ! reinitiialize distribution function do k=0,8 ! if(flag(i+cx(k),j+cy(k)).eq.3) then t1=aveuu*aveuu+aveuv*aveuv t2=aveuu*cx(k)+aveuv*cy(k) f(k,i,j)=w(k)*averh*(1.+3.0*t2+4.50*t2*t2-1.50*t1) enddo |
Харилцах үеийн нүдний хувьд хязгаарын нөхцлийг шийдвэрлэх хэрэгтэй. Үүнийг дараах байдлаар кодолж болно.
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 | ! boundary for interface ! normal unit vector (marchin square) do i=1,n-1 do j=1,m-1 if(flag(i,j).eq.2) then nx(i,j)=0.5*(e(i-1,j)-e(i+1,j)) ny(i,j)=0.5*(e(i,j-1)-e(i,j+1)) end if enddo enddo ! boundary for normal vector do i=0,n nx(i,0)=nx(i,1) nx(i,m)=nx(i,m-1) ny(i,0)=nx(i,1) ny(i,m)=nx(i,m-1) enddo do j=0,m nx(0,j)=nx(1,j) nx(n,j)=nx(n-1,j) ny(0,j)=nx(1,j) ny(n,j)=nx(n-1,j) enddo ! equilibrium function for air substances do i=0,n do j=0,m if(flag(i,j).eq.2) then t1=u(i,j)*u(i,j)+v(i,j)*v(i,j) do k=1,8 t2=u(i,j)*cx(k)+v(i,j)*cy(k) faq(k,i,j)=w(k)*1.0*(1.+3.0*t2+4.50*t2*t2-1.50*t1) enddo endif enddo enddo ! reconstruction for distribution function do i=0,n do j=0,m if(flag(i,j).eq.2) then do k=1,8 op=opp(k) if(cx(k)*nx(i,j)+cy(k)*ny(i,j).ge.0 & .or.flag(i+cx(k),j+cy(k)).eq.3) then f(op,i,j)=faq(k,i,j)+faq(op,i,j)-f(k,i,j) endif enddo endif enddo enddo print*, 'Boundary condition has performed' |
Тооцон бодогч: Усны барилгын инженер Б.Аюурзана.
No comments:
Post a Comment