Газрын доорхи усыг ЛБА-аар загварчилах үндэс /Groundwater flow simulation with LBM/

Газрын доорхи усыг ЛБА-аар загварчилах үндэс

Газрын доорхи усны урсгалыг энгийнээр /бүхнийг энгийнээр бодитоор гэдэг бол зөв бодох үндэс/ авч үзвэл өмнөө маш олон хөрсний үрлэн саадтай шингэний урсгал гэсэн үг. Урсгалын хурд маш бага байх ба гол хөдөлгөгч хүч нь хүндийн хүч, капиляр хүч байна. Газрын доорхи усны хөдөлгөөн бүх л шингэний урсгалыг илэрхийлэгч Навьер-Стокесийн тэгшитгэлд захирагдах ба үүнд эсэргүүцлийн илэрхийллийг оруулахад л хангалттай. Дарсийн хууль бол энэ эсэргүүцлийн нэг ойролцоолол юм. Сүвэрхэг орчны урсгал нь усны барилгад ихээхэн тохиолдоно. Шороон боомтын /аливаа шороон байгууламж/ дотуурх болон бууриар шүүрэх шүүрэлтийн урсгал бол энэ төрлийн барилгын тооцооны болоод аюулгүй ажиллагааны чухал хэсэг юм. Бетон байгууламжийн дотуур ч хүртэл шүүрэлт явагдаж байдаг. Бетоны ан цавд ус орж хөлдөн түүнийг эдвэх нь судлах, сэргийлэх ёстой асуудлуудын нэг юм. Уснаас мөс үүсэх процессийг бодчихвол монголд тулгамдаж байгаа нилээд хэдэн асуудлыг шийдчих боломж бүрдэнэ. Тухайлбал асфальтан замын хагарал, хөөлт, хотойлт гэх мэт. 

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

Доорхи илтгэлийн pdf хувилбарыг эндээс авна уу.

Тэкстэн дунд байгаа одууд бол томъёоны аль нэг гишүүдийн тайлбар гэдгийг анхааруулъя. Завгүй болхоор оруулж амжихгүй юм.

СҮВЭРХЭГ ОРЧИН ДОТУУРХ ШИНГЭНИЙ УРСГАЛЫГТООЦОН БОДОХ

Бадарчийн Аюурзана, Эрдэнэбаярын Тунгалагтамир

Япон улсын Нагаока Технологийн их сургууль

ХУРААНГУЙ

Шинжлэх ухааны ихэнхи салбар тэр дундаа инженерчлэлийн салбарын судалгааны чухал чиглэл бол тооцон бодох ухаан юм. Ус, шингэнтэй холбогдсон инженерчлэлийн салбарт уг ухаан нь тооцон бодох шингэний динамик хэмээн нэрлэгдэж салбарын чухал асуудлуудыг шийдэхэд үнэтэй хувь нэмрээ оруулж байна. Энэхүү илтгэлд харьцангуй шинэвтэр шингэний хөдөлгөөнийг тооцон бодох арга болох Латтис Больцманы аргыг (цаашид ЛБА гэх) сүвэрхэг орчин дундуурх шингэний урсгалд хэрхэн ашиглаж байгаа тухай дурьдах юм.  ЛБА нь бусад уламжлалт тооцон бодох аргуудтай (төгсгөлөг аргууд) харьцуулбал сүвэрхэг орчныг загварчлахад маш тохиромжтойн дээр, хэд хэдэн ялгаатай замаар тооцон бодох боломжийг бүрдүүлж өгч байгаа нь хэмжээсийн янз бүрийн түвшинд газрын доорхи усны хөдөлгөөн, түүнд ууссан бодисын тархалт/дисперси зэргийг загварчлах боломжтойг харуулж байна. Энэ илтгэлийн үр дүнд гарсан тооцооны тохирол болон өмнөх ижил төстэй судалгааны материалуудаас үзэхэд ЛБА нь сүвэрхэг орчин дотуурх урсгалын динамик шинж чанарыг загварчлах өндөр бүтээмжтэй, нарийвчлалтай арга гэдгийг баталж байна.

Түлхүүр үг: Тооцон бодох шингэний динамик, Латтис Больцманы арга, сүвэрхэг орчин дотуурх шингэний хөдөлгөөн.

ОРШИЛ

Тооцон бодох шингэний динамикийн тооцооны аргуудыг ерөнхийд нь базаж харвал хэлхээст болон чөлөөт аргууд гэж ангилагдана. Хэлхээст гэдэг ангилалд загварчлах объектийг геометрийн хувьд орон зайн хэлхээст хувааж түүний зангилаа, шилбэ гэх мэт дээр шингэний физик хэмжигдэхүүний хугацаан дахь өөрчлөлтийг тооцон боддог Эйлэрийн аргууд багтана. Чөлөөт аргад харьцангуй хөдлөх боломжтой цэгэн дээрх шингэний шинж чанаруудыг тооцон бодох Лагранжийн хандлагатай аргууд багтана. Газрын доорхи усны хөдөлгөөнийг хэлхээст аргаар загварчилсан маш олон хэрэглээний болон чөлөөт програм хангамжууд өргөн дэлгэр ашиглагдаж байна.  Бидний авч үзэх ЛБА нь мөн хэлхээст аргын ангилалд багтах боловч бусад хэлхээст арга болох төгсгөлөг ялгавар болон төгсгөлөг эзлэхүүний аргуудаас зарим нэг давуу талтай учир судлаачдын анхаарлыг маш их татаж байна. Дурьдваас, төвөгтэй геометрийн бүтэцтэй обьектыг хялбархан загварчлах, том хэмжээст загварчлалыг параллель програмчлалын зарчмаар богино цаг хугацаанд тооцох, боловсруулалтын хувьд маш энгийн, хэмжээсгүй орчинд тооцоо хийдэг гэх мэт болно. Энэ илтгэлд уг аргын зөвхөн сүвэрхэг орчны урсгалтай холбоотой тооцон бодох аргачлалыг, хоёр хэмжээст орчинд хэрхэн ашиглах тухай авч үзэх учир практик ашиглалтаас дурьдахгүй. Загварыг зөв боловсруулж (код бичих) зүгшрүүлсэн тохиолдолд практик хэрэглээнд нэвтрүүлж, сав газар эсхүл усны ордын хөдөлгөөнийг загварчилах нь дараагийн алхам бөгөөд төвөггүй гэж бодож байна.

АРГА, АРГАЧЛАЛ

1.    Латтис Больцманы аргын тухай

ЛБА-ын онолын үндэслэл нь идиал хийн кинетик дээр тогтдог юм [1]. Хийн кинетикийг илэрхийлэх Больцманы тэгшитгэлийг тодорхой чиглэлийн хурднуутай сүлжээн (lattice) дээр ойролцоолсон тэгшитгэлийг Сүлжээний Больцманы тэгшитгэл гэнэ [2].

Дээрх тэгшитгэлд буй f хэмжигдэхүүн нь хийн молекулын түгэлтийг илэрхийлэх ба тэгшитгэлийн баруун талын илэрхийлэл * нь молекулуудын мөргөлдөөнийг илэрхийлэх оператор юм. Гаднаас хүч * үйлчилж буй Больцманы тэгшитгэлийг сүлжээн дээр ойролцоолоход хамгийн бэрх гишүүн нь дээрх мөргөлдөөний оператор байх бөгөөд хамгийн хялбар ойролцооллыг дан тайвшралттай шугаман байхаар тодорхойлсон байна [3]. Энэ ойролцоололтой (ВГК гэж товчилж нэрлэнэ) сүлжээний тэгшитгэл нь дараах байдлаар илэрхийлэгдэнэ.

Үүнд: * нь * байрлалд байх * хугацаан дахь түгэлтийн функ, * нь * байрлалд байх * хугацаан дахь түгэлтийн функц, * нь түгэлтийн функцыг тэнвцэрийн төлөврүү чиглүүлэх тэнцвэрт түгэлтийн функц, * нь шингэнд үйлчилж буй гадаад хүч, * нь тайвшралтын хэмжээсгүй хугацааг илэрхийлнэ. Тэнцвэрийн нөхцөлрүү хөтлөх тайвшралтын хэмжээсгүй хугацаа нь хурдны оронд сүлжээний * зунгааралттай шууд хамааралтай байна. Түгэлтийн функц нь * хугацааны алхамтайгаар сүлжээний чиглэлийн * хурдны дагуу нүүж, шинэ сүлжээнд очихдоо бусад сүлжээнээс ирсэн түгэлтийн функцуудтайгаа мөргөлдөх үйл явцад орно. Энэхүү илтгэлд онолын тайлбарыг орхиж хэрэглээ болон боловсруулалтад түлхүү тайлбарыг өгч байгаа юм. Хоёр хэмжээст орчинд ихэвчлэн хэрэглэгддэг сүлжээ нь 8 чиглэлд эгэл хурдтай байх ба D2Q9 гэж нэрлэгддэг (Зураг 1).

Зураг 1. Зураг (а) нь D2Q9 сүлжээний байгуулалт, сүлжээний тухайлсан хурднуудын чиглэл, тэмдэглэгээ, Зураг (б) нь хуучин хугацаан дахь түгэлтийн функцуудыг илэрхийлэх ба зураг (в) нь * хугацааны алхамын дараах тэдгээрийн нүүлтийг, зураг (г) нь шинэ байрлал дээр мөргөлдөөний операторыг тооцсоны дараах шинэ түгэлтийн функцуудыг илэрхийлнэ.

Сүлжээний тухайлсан хурднууд нь жигд сүлжээний хувьд утга нь ямагт нэгж байна. Тэнцвэрт түгэлтийн функц болон тайвшралтын хугацаа нь тооцооны чухал үүргийг гүйцэтгэх ба загварчилж буй шингэний хөдөлгөөний онцлог шинж чанарыг агуулсан байх учиртай. Сүлжээний Больцманы тэгшитгэлд Чапман-Энскогийн задаргааг хийснээр нил орчны шингэний механик дахь (макро хэмжээст) Навьер-Стокесийн тэгшитгэлийг сэргээж болох ба энд тохирох тэнцвэрт түгэлтийн функц нь дараах байдлаар тодорхойлогдоно.

Үүнд * нь сүлжээний жинлэх итгэлцүүр, * нь сүлжээний дууны хурд, * болон * нь макро хэмжээст дахь шингэний нягт болон хурднууд юм. Эдгээрийг түгэлтийн функцын тэг ба нэгдүгээр эрэмбийн моментуудын нийлбэрээр тодорхойно.

2.      Сүвэрхэг орчин дундуурх шингэний хөдөлгөөний Латтис Больцманы загвар

Сүвэрхэг орчны хэмжээ, хоосон орон зайн нөлөөлөл, ач холбогдлыг авч үзсэн хоёр янзын аргачлал (зураг 2) ЛБА-д хэрэглэгдэж байна.

Зураг 2. Сүвэрхэг орчин дахь урсгалыг загварчлах ЛБ-ны аргачлал. Эхний зураг нь сүвний хэмжээст аргачлалыг харуулах бол удаах нь эзлэхүүний элементар төлөөлөлт аргачлалыг харуулж байна.

Хэрэв сүвэрхэг орчны бүтэц, нүх сүвийн хэмжээ, түүний геометрийн онцлог нь маш чухал, шингэний урсгалыг микро орчинд судлаж үзэх шаардлагатай бол сүвний хэмжээст (pore-scale approach) аргачлалыг авч үзнэ. Энэ аргачлал нь дээр тайлбарласан ердийн ЛБА-тай яг ижилхэн бөгөөд хязгаарын нөхцөлд маш нямбай хандах шаардлага гардаг. Хэрэв сүвний хэмжээст тооцоолол шаардлагагүй, том орон зайд тооцоолох шаардлагатай бол эзлэхүүний элементар төлөөллийн (Representative elementary volume approach) аргачлалыг хэрэглэнэ.  Альч аргачлал нь хөрсний болоод газрын доорхи усны хөдөлгөөний загварчлалд чухал юм.

3.      Нүх сүвний хэмжээст аргачлал

Нүх сүвний хэмжээ нь хэдэн нано метрийн хэмжээнд яригдах боловч ЛБА-ыг хэрэглэхэд энгийн аргачлалтай яг адил хэрэглэнэ. ЛБА нь Сүлжээний Больцманы тэгшитгэлийг хоёр алхамаар тооцно.

Нүүх алхамаар тодорхойлогдсон түгэлтийн функц нь тэгшитгэл 6 – ийн баруун гар талд орлуулагдаж шинэ түгэлтийн функц * – ийг өгнө.

Сүвэрхэг орчны шингэнтэй харьцах гадаргууг загварчилахдаа тохирох хязгаарын нөхцлийг өгнө. Ихэнхи хатуу гадаргуу нь барзгарын нөхцлөөс хамаарч гулсах эсхүл үл гулсах хязгаарын нөхцлийг авна. ЛБА-д хатуу гадаргуугийн хязгаарын нөхцлийг түгэлтийн функцын буцаж ойх схемээр шийддэг. Түгэлтийн функц нь хатуу гадаргуугийн зангилаанд нүүж ирсэн бол ирсэн зүгрүүгээ буцаж ойх нь үл гулсах хязгаарын нөхцлийг илэрхийлнэ. Сүвэрхэг орчны хана геометрийн хувьд бэрх байх боловч ЛБА-д буцаж ойх схемийг хэрэглэсэнээр маш хялбархан шийдэгддэг. Сүвэрхэг орчнийг зөв зурагласан гурван хэмжээст эсвэл хоёр хэмжээст геометрийн өгөгдөлтэй үед энэ аргачлалаар хялбархан тооцож болохыг үр дүн хэсгээс үзээрэй.

4.      Эзлэхүүний элементар төлөөллийн аргачлал

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

а. Сүвэрхэгшилтэнд суурилсан арга

Нэгэнт орчны нүх сүвийг эзлэхүүний түвшинд орлуулах төлөөлөгчөөр сүвэрхэгшилтийг сонгосон бол түгэлтийн функцууд нь эзлэхүүнээр дундажлагдах ёстой. Сүвэрхэг орчны шугаман эсэргүүцэл болох Дарси, шугаман бус эсэргүүцэл болох Форчшэймэрийн тэгшитгэлийг хангасан макро хэмжээс дахь урсгалын тэгшитгэлийг эзлэхүүнээр дундачилсан Больцманы тэгшитгэлээс сэргээвэл [4]:

Тэгшитгэл 7 нь урсгал тасралтгүйн тэгшитгэлийг, тэгшитгэл 8 нь хөдөлгөөний тоо хэмжээний өөрчлөлтийн тэгшитгэлийг сүвэрхэг орчинд илэрхийлнэ. Гадаад хүчинд  уг орчний шугаман ба шугаман бус эсэргүүцлийг хамтад нь багтаана.

Үүнд * нь орчны сүвэрхэгшилт, * – шүүрэлтийн итгэлцүүр, * нь сүвэрхэг орчны геометр бүтцийн функц [5],  * биеийн болон гадаад хүчийг илэрхийлнэ.  Тэгшитгэл 7 ба 8-аар тодорхойлогдох сүвэрхэг орчны урсгалд дараах хэмжээсгүй тоонууд (Рейнольдийн тоо, Дарсийн тоо, Зунгааралтын харьцаа) тодорхойлогдох ба ЛБА-д шилжүүлэгчид болно.

Рейнольдсын тоо болон Дарсийн тоо хэт бага байхад эсэргүүцлийн шугаман бус илэрхийллийн нөлөөлөл бага байх учир зарим тохиолдолд орхигддог. Эзлэхүүнээр дундажлагдсан урсгалыг загварчлахад сүлжээний Больцманы тэгшитгэл нь стандарт тэгшитгэлтэй (тэгшитгэл 2) ижил байна. Харин хүчний илэрхийлэл, тэнцвэрт түгэлтийн функцууд нь дараах байдлаар тооцогдоно.

Зарим схемд хүчийг тооцохдоо хурдны тооцоонд нэмэх ба энэ тохиолдолд тэгшитгэл 4 дахь хурдны томъёо дараах байдалтай болно.

Хүчний илэрхийллийг тэгшитгэл 11-ээр тооцоход хурдыг мэдсэн байх шаардлагатай ба түр хугацааны хурдыг тодорхойлох аргачлалыг лавлах [6]–аас үзэж болно.

б. Хатуу хэсгийн индекст суурилсан арга

Шингэний урсгалын момент нь хатуу хэсгүүдийн эсэргүүцлийн улмаас алдагдах ба энэ санаанд суурилж Бринкманы өргөтгөсөн Дарсийн тэгшитгэлийг хангах ЛБА-ыг боловсруулсан байна [7]. Уг аргад сүвэрхэгшилээс хамаарах хатуу хэсгийн индекс * хэмээх хэмжигдэхүүнийг ашиглах ба уг хэмжигдэхүүн нь шүүрэлтийн итгэлцүүртэй * гэж холбогдоно. Урсгалын дараах түгэлтийн функцыг (тэгшитгэл 5-аас) * гэж тэмдэглэвэл мөргөлдөөн (тэгшитгэл 6) нь дараах байдлаар бичигдэнэ.

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

Үүнд: * нь * – ын эсрэг чиглэлийг үзүүлэх ба * бол урсгалд моментийн өөрчлөлт байхгүй буюу чөлөөт урсгалтай, * үед хатуу хэсгийн нөлөөлөл тооцогдож урсгал моментоо алдах болно. Зарим судлаачид үүнийг магадлалт ойлтын арга гэж нэрлэдэг. Дээрх тайлбарласан гурван аргаар гидравликийн зарим жишиг бодлогыг бодсон үр дүнгээс танилцуулъя.

ТООЦОН БОДОЛТЫН ҮР ДҮН

1.      Нүх сүвний хэмжээст ЛБА-ын үр дүн

Бодит хөрсний бүтцийг үзүүлэх 2 хэмжээст геометрийн өгөгдөл хараахан байхгүй учир санамсаргүй түгэлтээр байгуулагдсан сүвэрхэг орчинг авч үзье. Тооцооны параметрүүд болон анхны нөхцлүүд бүгд хэмжээсгүй байдлаар өгөгдөх учиртай ба бодит бодлого, ЛБ загварчлал хоорондоо хэмжээсгүй тоогоор холбогдох ёстой.

Зураг 3. Стандарт ЛБА-аар тооцсон нүх сүвний хэмжээст загварчлалын үр дүн. Эхний зураг нь санамсаргүй түгэлтээр бий болсон сүвэрхэг орчны геометр бүтцийг харуулна. Удаах зураг нь уг сүвэрхэг орчин дундуур урсгах шингэний урсгалын хурдны хэвтээ бүрдүүлэгчийг хурдны вектор оронтой хамтатган харуулж байна.

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

2.      Сүвэрхэгшилтэнд суурилсан арга

Боловсруулалтыг зөв хийсэн эсэх, тооцооллын нарийвчлалыг баталгаажуулах хамгийн эхний алхам бол гидравликийн жишиг бодлогын аналитик, эсхүл ялгаатай аргачлалаар бодсон үр дүнтэй харьцуулах юм. Сүвэрхэгшилтэд суурилсан аргыг  сүвэрхэгшилттэй материалаар дүүргэгдсэн хоёр хавтгай хана дотуур урсах усргал болон хөндийн урсгалыг загварчилж төгсгөлөг элементийн аргын үр дүнтэй [8] харьцуулав (зураг 4). Пуазейлийн урсгалын хувьд манай тооцооны үр дүн хурдны муруйн төв хэсэгтээ лавлахтай [5] тохирч байгаа боловч хязгаарын орчимд том зөрүүтэй байна. Харин хөндийн урсгалын хувьд маш ойролцоо шийдлийг өгч байна. Үр дүнгийн зөрүү нь хязгаарын нөхцлийн өгөгдлөөс хамаарч байгааг тодруулсан бөгөөд сүлжээний параметрийн зохистой харьцаа мөн чухал байлаа. ЛБА-аас гарч байгаа үр дүнг тодорхой масштабын хуулиар бодит хэмжээст үр дүнлүү шилжүүлдэг.

Зураг 4. Сүвэрхэг орчинтой үе дотуурх хүчилсэн Пуазейлийн урсгал. Дарсийн тоо тогтмол, Рейнольдсын өөр тоотой үеийн үр дүнгүүдийн харьцуулалт. Баруун талын зураг нь босоо хурдны оронг харуулна. Материалын сүвэрхэгшилт * байна.

Зураг 5. Хөндий нь * сүвэрхэгшилтэй, Дарсийн тоо тогтмол, Рейнольдсын өөр тоонуудтай байх үеийн үр дүнг төгсгөлөг ялгаварын аргатай харьцуулсан график.

3.      Хатуу хэсгийн индекст суурилсан арга

Санамсаргүй түгэлтээр байгуулагдсан сүвэрхэг орчинг ашиглаж 100х50 хэмжээтэй сүлжээн дээр сүвэрхэг орчныг эл аргаар мөн тооцож үзэв. Хэрэв яг ижил геометрийг стандарт болон магадлалт ойлтын аргаар бодож харьцуулбал маш ойролцоо үр дүнг үзүүлэх боловч стандарт ЛБА нь хатуу хэсэг тус бүрд буцаж ойх дүрмийг хэрэгжүүлдэг учир урсгалын эсэргүүцэл ихтэй үр дүнг өгнө.

Зураг 6. Рейнольдсын тоо Re=2 үеийн сүвэрхэг орчин дотуурх шингэний урсгал.

ДҮГНЭЛТ ХЭЛЭЛЦҮҮЛЭГ

ЛБА аргаар сүвэрхэг орчныг загварчлах аргуудыг тайлбарлаж, үндсэн тэгшитгэл, томъёоллыг тавьж, зарим нэг жишиг бодлогийг бодож төгсгөлөг ялгаварын аргаар бодсон үр дүнтэй харьцуулж тооцооны нарийвчлалыг шалгаж баталгаажуулав. Илтгэлд дурьдагдсан сүвэрхэг орчны урсгалыг загварчлах ЛБА нь энгийн бөгөөд маш үр ашигтайгаар шингэний урсгалыг микро ба макро хэмжээнд загварчлах учир сайтар судалж хэрэглээнд нэвтрүүлэх шаардлагатай. Дан ганц усны хөдөлгөөнөөс гадна усанд ууссан бодисын тархалт/дисперсийг мөн загварчлах, термонамиктай холбогдох бүрэн боломжтой. Газрын доорхи уснаас гадна шингэн ашигт малтмалын судалгаанд мөн өргөн ашиглагдаж байна. Энэ арга нь бусад тооцооллын аргатай харьцуулахад боловсруулах тэгшитгэл математикийн хувьд хүндрэлгүй, үйлчилэх хүрээ өргөнтэй учир судлаачдын анхаарлыг ихээр татаж байна. Дурьдагдсан тооцооны бодлогууд дахь параметр, анхны нөхцөл ба хязгаарын нөхцөл зэргийг багтаагүйн улмаас дурьдаагүй боловч лавлаж асуувал нарийн мэдээлэл олж авах боломжтой.

Ашигласан ном, сэтгүүл

[1]	S. Chen and G. D.Doolen, "Lattice Boltzmann Method for Fluid Flows," Annual Review of fluid mechanics, vol. 30, pp. 329-364, 1998.
[2]	A. A. Mohamad, Lattice Boltzmann method: fundamentals and engineering applications with computer codes, Springer Science and Bussiness Media, 2011.
[3]	P. L. Bhatnagar, E. P.Gross and M. Krook, "A model for collision processes in gases. I. Small amplitude processes in charged and neutral one component system," Physical review, vol. 94, no. 3, p. 511, 1954.
[4]	P. Nithiarasu, K. Seetharamu and T. Sundararajan, "Natural convective heat transfer in a fluid saturated variable porosity medium," International Journal of Heat and Mass Transfer, vol. 40, no. 16, pp. 3955-3967, 1997.
[5]	Z. Guo and T.S.Zhoa, "Lattice Boltzmann model for incompressible flows through porous media," Physical Review E, vol. 66, p. 036304, 2002.
[6]	Z. Guo and C. Shu, Lattice Boltzmann Method and its Application in engineering, Singapore: World scientific Publishing, 2013.
[7]	O. Dardis and J. McCloskey, "LAttice boltzmann scheme with real numbered solid density for the simulation of flow in porous media," Physical review E, vol. 57, no. 4, pp. 4834-4837, 1998.
[8]	U.Ghia, K.N.Ghia and C.T.Shin, "High-Re Solutions for Incompressible flow using Navier-Stokes Equations and a Multigrid Method," Journal of computational physics, vol. 48, pp. 387-411, 1982.