Bir muncha vaqt oldin men OpenFOAM raqamli modellashtirish kutubxonasida sodir bo'ladigan operatsiyalar va jarayonlarning tavsifini qidirayotgan edim. Men cheklangan hajm usulining ishlashining ko'plab mavhum tavsiflarini, klassik farq sxemalarini va turli xil fizik tenglamalarni topdim. Men batafsilroq bilmoqchi edim - bu qiymatlar falon va falon chiqish faylida falon iteratsiyada qaerdan paydo bo'lgan, fvSchemes, fvSolution sozlamalari fayllaridagi ma'lum parametrlar ortida qanday iboralar bor?
Bunga qiziqqanlar uchun - ushbu maqola. OpenFOAM yoki unda amalga oshirilgan usullar bilan yaxshi tanish bo'lganlar - shaxsiy xabarda topilgan xatolar va noaniqliklar haqida yozing.

Habré-da OpenFOAM haqida bir nechta maqolalar allaqachon mavjud edi:

Shuning uchun men bu "sonli simulyatsiya uchun ochiq (GPL) platformasi, cheklangan hajm usuli yordamida qisman differensial tenglamalarni echish bilan bog'liq simulyatsiyalar uchun mo'ljallangan va uzluksiz mexanikadagi muammolarni hal qilish uchun keng qo'llaniladigan" ekanligiga to'xtalib o'tirmayman.

Bugun men OpenFOAM-da hisob-kitoblar paytida sodir bo'ladigan operatsiyalarni tasvirlash uchun oddiy misoldan foydalanaman.

Shunday qilib, geometriyani hisobga olgan holda - tomoni 1 metr bo'lgan kub:

Bizning oldimizda ma'lum bir skalyar maydonning (harorat, moddaning miqdori) oqimining tarqalishini modellashtirish vazifasi turibdi, bu jism hajmi ichidagi quyidagi transport tenglamasi (1) bilan beriladi.

(1)
,

Skayar miqdor, masalan, haroratni [K] yoki ma'lum bir moddaning konsentratsiyasini ifodalasa va moddaning o'tkazilishini ifodalaydi, massa oqimi [kg/s].

Bu tenglama, masalan, issiqlik tarqalishini modellashtirish uchun ishlatiladi
,
Bu erda k - issiqlik o'tkazuvchanligi va harorat [K].

Divergentsiya operatori aslida

operator.
Eslatib o'taman, nabla operatori (Gamilton operatori) mavjud bo'lib, u quyidagicha yozilgan:
,

Bu yerda i, j, k - birlik vektorlar.
Agar nabla operatorini vektor kattalikka skalyar ko'paytirsak, bu vektorning divergensiyasini olamiz:

"Fizika nuqtai nazaridan, vektor maydonining divergensiyasi kosmosdagi ma'lum bir nuqta bu maydonning manbasi yoki cho'kishi darajasining ko'rsatkichidir"

Agar siz nabla operatorini skalerga ko'paytirsangiz, siz ushbu skalerning gradientini olasiz:

Gradient skalyar kattalikning qaysidir yo'nalishda o'sishi yoki kamayishini ko'rsatadi.


Muammoning chegara shartlari quyidagicha: kirish yuzi, chiqish yuzi, qolgan yuzlari esa silliq devorlardir.

Kub hajmini cheklangan hajmlarga bo'lish

Bizning panjaramiz juda oddiy bo'ladi - biz kubni Z o'qi bo'ylab 5 ta teng hujayraga ajratamiz.

Ko'p formulalar

Cheklangan hajm usuli (1) integral shakldagi (2) har bir chekli hajm uchun qanoatlantirilishini ta'minlaydi.

(2)
,

Yakuniy hajmning geometrik markazi qayerda.

Yakuniy hajm markazi


(2) ifodaning birinchi hadini quyidagicha soddalashtiramiz va o'zgartiramiz:

(2.1) (HJ-3.12)*

Ko'rib turganingizdek, biz skalyar miqdor cheklangan hajm ichida chiziqli ravishda o'zgaradi deb taxmin qildik va chekli hajm ichidagi bir nuqtadagi miqdorning qiymatini quyidagicha hisoblash mumkin:

Ifodaning ikkinchi hadini (2) soddalashtirish uchun umumlashtirilgan Gauss-Ostrogradskiy teoremasidan foydalanamiz: vektor maydonining hajm bo'yicha divergensiyasining integrali berilgan hajmni chegaralovchi sirt orqali vektor oqimiga teng. Inson tilida "cheklangan hajmga kiruvchi/chiqadigan barcha oqimlarning yig'indisi ushbu cheklangan hajmning yuzlari orqali o'tadigan oqimlarning yig'indisiga teng":

(2.3)
,

Ovozni cheklaydigan yopiq sirt qayerda,
- hajmdan normal bo'ylab yo'naltirilgan vektor.

Vektor S



Cheklangan hajm tekis yuzlar to'plami bilan cheklanganligini hisobga olsak, (2.3) ifodani sirt ustidagi integrallar yig'indisiga aylantirish mumkin:

(2.4) (HJ-3.13)
,

Bu erda yuzning markazida o'zgaruvchining qiymatini ifodalaydi,
- maydon vektori, yuzning markazidan chiqadigan, hujayradan uzoqqa (mahalliy), pastroq indeksli hujayradan yuqori indeksli (global) katakchaga yo'naltirilgan.

S vektori haqida bir oz ko'proq

Xuddi shu vektor parametrlarini ikki marta saqlamaslik uchun, chunki Ko'rinib turibdiki, ikkita qo'shni hujayralar uchun hujayra o'rtasidan uzoqqa yo'naltirilgan hujayralar orasidagi chekka normal vektor faqat yo'nalish belgisi bilan farq qiladi. Shuning uchun chekka va hujayra o'rtasida ega-qo'shni munosabatlari yaratilgan. Agar maydon vektori (global, indeksi pastroq katakdan kattaroq indeksli katakgacha bo‘lgan musbat yo‘nalish) hujayra markazidan FROM ni ko‘rsatsa, yacheyka va vektor o‘rtasida, aniqrog‘i katak va hujayra o‘rtasida shunday munosabat bor. yuz, egani bildiradi). Agar ushbu vektor ko'rib chiqilayotgan hujayra ichiga ishora qilsa, u holda qo'shni. Yo'nalish qiymat belgisiga ta'sir qiladi (egasi uchun + va qo'shni uchun) va bu yig'ishda muhim, pastga qarang.

Farq sxemalari haqida

Yuzning markazidagi qiymat qo'shni hujayralar markazlaridagi qiymatlar orqali hisoblanadi - bu ifodalash usuli farq sxemasi deb ataladi. OpenFOAM-da farq sxemasining turi faylda ko'rsatilgan /system/fvSchemes:

DivSchemes (standart yo'q; div(phi, psi) Gauss chiziqli; )

Gauss- markaziy farq sxemasi tanlanganligini bildiradi;
chiziqli- hujayralar markazlaridan yuzlar markazlarigacha bo'lgan interpolyatsiya chiziqli ravishda sodir bo'lishini anglatadi.

Faraz qilaylik, bizning skalyar miqdorimiz chekli hajm ichida markazdan qirralarga chiziqli ravishda o'zgaradi. Keyin yuzning markazida taxminiy qiymat quyidagi formula bo'yicha hisoblanadi:

Og'irliklar qayerda va sifatida hisoblanadi

Hujayra hajmlari qayerda.
Egri hujayralar holatlari uchun, taxminiy og'irliklarni hisoblash uchun murakkabroq formulalar mavjud.

Shunday qilib, hujayraning chekka markazlaridagi phi_f qiymatlari hujayra markazlaridagi qiymatlar asosida hisoblanadi. Gradient qiymatlari grad(phi) phi_f qiymatlari asosida hisoblanadi.
Va bu butun algoritm quyidagi psevdokod shaklida taqdim etilishi mumkin.
1. Biz chekli hajmli gradientlar massivini e'lon qilamiz, uni nollar bilan ishga tushiramiz 2. Biz barcha ichki yuzlardan o'tamiz (ular chegara bo'lmagan) > Biz flux_f = phi_f*S_f ni hisoblaymiz. Hujayra sentlarida phi qiymatlari asosida phi_f qiymatlarini hisoblang > Ega elementining gradientiga flux_f va qo'shni elementning gradientiga -flux_f qo'shing 3. Barcha chegara yuzlari bo'ylab takrorlang > Flux_f = phi_f*S_f ni hisoblang > Ega elementining gradientiga flux_f qo'shing (qo'shni -chegara yuzlari elementga ega emas) 4. Keling, barcha elementlarni ko'rib chiqamiz > Olingan gradient yig'indisini element hajmiga bo'ling.

Vaqt namunasi

(2.1) va (2.4) ni hisobga olgan holda (2) ifoda quyidagi shaklni oladi:

(3)

Cheklangan hajm usuliga ko'ra, vaqtni diskretlashtirish amalga oshiriladi va (3) ifoda quyidagicha yoziladi:

(4)

Keling, integrallashamiz (4):

(4.1)

Keling, chap va o'ng tomonlarni ajratamiz:

(5)

Namuna olish matritsasi uchun ma'lumotlar

Endi har bir chekli hajm uchun chiziqli tenglamalar tizimini olishimiz mumkin.

Quyida biz foydalanadigan panjara tugunlarining raqamlanishi keltirilgan.

Tugun koordinatalari /constant/polyMesh/punktlarida saqlanadi

24 ((0 0 0) (1 0 0) (0 1 0) (1 1 0) (0 0 0.2) (1 0 0.2) (0 1 0.2) (1 1 0.2) (0 0 0.4) (1 0 0.4) (0 1 0.4) (1 1 0.4) (0 0 0.6) (1 0 0.6) (0 1 0.6) (1 1 0.6) (0 0 0.8) (1 0 0.8) (0 1 0.8) (1 1 0.8) (0 0 1) (1 0 1) (0 1 1) (1 1 1))

Hujayralarning tugunlari-markazlarini raqamlash (50, 51 - chegara yuzlari markazlari):

Yuzning markaziy tugunlarini raqamlash:

Elementlar hajmi:

Hujayra yuzidagi qiymatlarni hisoblash uchun zarur bo'lgan interpolyatsiya koeffitsientlari. "e" pastki belgisi "hujayraning o'ng chetini" bildiradi. Ko'rinishga nisbatan to'g'ri, "Yacheykalar tugunlari-markazlarini raqamlash" rasmidagi kabi:

Namuna olish matritsasini shakllantirish

P = 0 uchun.
Miqdorning harakatini tavsiflovchi ifoda (5).

Har bir shakldagi chiziqli algebraik tenglamalar tizimiga aylantiriladi:

Yoki yuzlardagi nuqtalar indekslariga ko'ra

Hujayraga/xujayradan barcha oqimlarni yig'indi sifatida ifodalash mumkin

Bu erda, masalan, E katakning markaziy nuqtasida oqimning linearizatsiya koeffitsienti,
- yuzning markaziy nuqtasida oqimning linearizatsiya koeffitsienti,
- chiziqli bo'lmagan qism (masalan, doimiy).

Yuzlarning raqamlanishiga ko'ra, ifoda quyidagi shaklni oladi:

P_0 elementi uchun chegara shartlarini hisobga olgan holda, chiziqli algebraik tenglamani quyidagicha ifodalash mumkin.

...oldin olingan koeffitsientlarni almashtiring...

Kirish "a" dan oqim hujayra ichiga yo'naltiriladi va shuning uchun salbiy belgiga ega.

Bizning nazorat ifodamizda diffuziya atamasi bilan bir qatorda vaqt muddati ham mavjud, ammo yakuniy tenglama shunday ko'rinadi.

P = 1 uchun.

P = 4 uchun.

Chiziqli algebraik tenglamalar tizimi (SLAE) matritsa shaklida ifodalanishi mumkin

A(i,j) === 40,5 0,5 0 0 0 -0,5 40 0,5 0 0 0 -0,5 40 0,5 0 0 0 -0,5 40 0,5 0 0 0 -0,5 40,5

Psi = o'lchamlar; ichki maydonning bir xil bo'lmagan ro'yxati 5(0,0246875 0,000308546 3,85622e-06 4,81954e-08 5,95005e-10);

Buning asosida vektor uchun qiymatlar olinadi

Keyin vektor SLAE ga almashtiriladi va vektorni hisoblashning yangi iteratsiyasi sodir bo'ladi.

Va shunga o'xshash kelishmovchilik kerakli chegaralarga yetguncha davom etadi.

Havolalar

* Ushbu maqoladagi ba'zi tenglamalar Jasak Hrvoje (HJ - tenglama raqami) dissertatsiyasidan olingan va agar kimdir ular haqida ko'proq o'qishni xohlasa (

Ilgari, bir qator raqamli usullar uchun boshlang'ich nuqta bo'lib xizmat qilgan subdomain usuli eslatib o'tilgan. Bunday usullardan biri chekli hajm usulidir. Xuddi shu usul boshqa keng tarqalgan sinf - integral usullarning vakili. Subdomain usulini belgilashning klassik shaklidan hisoblash sohasining subdomenlarga bo'linishi va subdomen ustidagi qoldiqning integratsiyasi olinadi. Farqi, taxminiy (sinov) funktsiyasining aniq qayd etilishining yo'qligi. Ammo, avvalgidek, biz har bir subdomaindagi tenglamani "aniq" hal qilishga harakat qilamiz. Shuning uchun, asl tenglama subdomen ustida birlashtirilgan. Integral usullar shundan iboratki, avval differensial tenglamaning integrali olinadi va tenglamani yozishning integral shakli olinadi. Keyinchalik bu shakldagi tenglama alohida panjara hujayralariga qo'llaniladi. Bunday holda, hujayralar va pastki hududlar bir xil narsadir.

Darhaqiqat, tenglamalarni yozishning integral shakli (fizika nuqtai nazaridan) differentsialdan ko'ra kengroq qo'llanish doirasiga ega. Gap shundaki, funktsiya uzilishlari mavjud bo'lganda, differentsial tenglamalar qo'llanilmaydi va ularning integral analoglari ishlashda, ishlashda va ishlashda davom etadi .... Afsuski, ular raqamli ravishda amalga oshirilganda, bu afzallik ba'zan yo'qoladi.

Qoida tariqasida, tenglamalardan integrallar oddiy va tushunarli jismoniy ma'noga ega. Masalan, uzluksizlik tenglamasini ko'rib chiqing. Asl differensial tenglama yoziladi

Uni S sirtga ega bo'lgan V hajm bo'yicha va vaqt o'tishi bilan t 0 dan t 1 gacha bo'lgan oraliqda integrallaymiz. Hosilalarni integrallashda biz Stokes formulasidan foydalanamiz (uning maxsus holatlari Green va Ostrogradskiy-Gauss formulalari deb ataladi). Natijada biz olamiz

Bu belgida birinchi ikkita integral o'rtasidagi farq ko'rib chiqilayotgan vaqt oralig'ida berilgan hajmdagi massaning o'zgarishini anglatadi. Va qo'sh integral ma'lum vaqt oralig'ida uni chegaralovchi sirt orqali ma'lum hajmga oqib o'tadigan massani ko'rsatadi. Tabiiyki, biz raqamli usullar haqida gapirganimiz sababli, bu integrallar taxminan hisoblanadi. Va bu erda chekli farqlar usulida ko'rib chiqilganlarga o'xshash yaqinlashish savollari boshlanadi.



Keling, eng oddiy holatlardan birini ko'rib chiqaylik - ikki o'lchovli to'rtburchaklar bir xil panjara. Cheklangan hajm usulida funktsiyalarning qiymatlari odatda panjara tugunlarida emas, balki hujayralar markazlarida aniqlanadi. Shunga ko'ra, indekslangan har bir yo'nalishdagi panjara chiziqlari emas, balki hujayralar qatlamlari (rasmga qarang).

j-1
j
j+1
k-1
k
k+1
A
B
C
D

Bu holda tenglamaning integral shakli quyidagicha yoziladi

Ko'rib turganingizdek, bu holda biz oddiy tenglamani oldik, uni chekli farqlar usuli yordamida ham yozishimiz mumkin. Bu barqarorlikni o'rganishning bir xil usullarini unga nisbatan qo'llash mumkinligini anglatadi. (Tezkor savol: bu sxema barqarormi?)

Ammo bizda bir xil narsa bo'lsa, unda butun bog'ni qurishga arziydimi? Eng oddiy hollarda, biz, albatta, hech qanday foyda olmaymiz. Ammo murakkabroq vaziyatlarda foydalar paydo bo'ladi. Birinchidan, yuqorida ta'kidlab o'tilganidek, bunday usullar (hatto shunday oddiy amalga oshirishda ham) uzilishlar va yuqori gradientli maydonlarni ancha yaxshi tasvirlaydi. Shu bilan birga, massa, impuls va energiyaning saqlanish qonunlarining bajarilishi kafolatlanadi, chunki ular har bir hujayrada kuzatiladi. Ikkinchidan, bu usullar tarmoqdagi turli xil suiiste'mollarga bardosh bera oladi. Hatto egri, notekis va tartibsiz panjaralar ham bu usullarni izdan chiqarib yubormaydi. Bu imtiyozlar, ayniqsa, chegara shartlari aniqlanganda seziladi.

j-1
j
j+1
k-1
k
k+1
A
B
C
D
E

Masalan, rasmda ko'rsatilgan holat uchun tenglamaning integral shakli shaklga ega bo'ladi

ya'ni, biz to'liq katak maydoni bo'yicha integralni olgan joyimizda, endi biz uni "kesilgan" maydondan o'tkazamiz, bu erda biz integralni to'liq chetidan oldik, endi biz uni uning qolgan qismidan olamiz. . Chegara qismidagi integral qo'shildi. Ammo chegara shartlaridan osongina topiladi. Xususan, agar devor orqali hech qanday massa oqimi ta'minlanmasa (shuningdek, sirtdan hech qanday massa olib tashlanmasa va / yoki biz devorda zaryadni yo'qotadigan ionlarning massa oqimini e'tiborsiz qoldiramiz), unda bunday integral oddiygina nolga teng bo'ladi. Energiya tenglamasining shunga o'xshash shaklida, qoida tariqasida, devor orqali oqimni hisobga olish kerak. Lekin chegara shartlaridan (agar ular to'g'ri o'rnatilgan bo'lsa) topish ham qiyin emas.

Buni mustahkamlash uchun impulsni saqlash tenglamalaridan biriga chekli hajm usulini qo‘llash qanday ko‘rinishini tasvirlab beraylik. Yagona zaryadlangan ionlar uchun tekis statsionar holatni olaylik. Biz yopishqoqlik va elastik to'qnashuvlarni e'tiborsiz qoldiramiz. Biz tenglamani olamiz

To'rtburchaklar to'r uchun (yuqoridagi rasmga qarang) biz olamiz

Bunday tenglamaning eng oddiy yaqinlashuvini quyidagicha yozish mumkin:

qisqartirishdan keyin formulani olamiz

Foydalanish chekli (nazorat) hajm usuli Keling, ikki o'lchovli statsionar issiqlik tenglamasi misolida ko'rsatamiz:

Guruch. 13. (31) tenglamani yechishda foydalaniladigan hisoblash panjarasi.

cheklangan hajm usuli

O'rtacha qiymat teoremasidan foydalanib, biz yozishimiz mumkin

,

Bu yerda Dx, DU – hujayra yuzlarining uzunliklari, x W – A katakchaning chap (“g‘arbiy”) chegarasining abssissasi, x E – o‘ng (“sharqiy”) chegarasining abssissasi, y N – ordinatasi. yuqori ("shimoliy") chegaraning, y S - pastki ("janubiy") chegaraning ordinatasi, S * - hujayraning o'rtacha issiqlik chiqarish tezligi. (32) chap tomonidagi lotinlar bo'yicha indeks (*), chegaralarning har birida issiqlik oqimlarini to'g'ri ko'rsatadigan tarzda aniqlangan o'rtacha qiymatlar sifatida ko'rib chiqilishi kerakligini ko'rsatadi. Ushbu holatni hisobga olgan holda, (32) ning diskret analogini qiyinchiliksiz olish mumkin [Patankar].

Shunday qilib, (32) tenglama A hujayra ichidagi issiqlik balansini (energiya saqlanish qonuni) tavsiflaydi. Hujayralar orasidagi issiqlik oqimlari to'g'ri tasvirlangan taqdirda, har bir nazorat hajmiga qo'llaniladigan (32) shakldagi tenglamalardan tashkil topgan tizim to'g'ri bo'ladi. butun hisoblash sohasi bo'ylab issiqlik balansini tavsiflang.

Paragrafning oxirida shuni ta'kidlash kerakki, alohida hollarda yuqorida tavsiflangan usullar bilan olingan hisoblash formulalari mos kelishi mumkin va egri chiziqli ortogonal bo'lmagan hisoblash tarmoqlaridan foydalanganda eng muhim farqlar paydo bo'ladi.

5. Diskret sxemalarning xossalari

5.1 Aniqlik

Aniqlik amaliy foydalanish uchun raqamli sxemaning maqbulligini tavsiflaydi. Diskret sxemaning to'g'riligini baholash juda qiyin vazifa bo'lib tuyuladi, chunki kontaktlarning zanglashiga olib keladigan xususiyatlari natijasida yuzaga keladigan xatolarni boshqa omillar natijasida yuzaga keladigan xatolardan ajratish deyarli mumkin emas (masalan, yaxlitlash xatolari, chegara va boshlang'ich shartlarni belgilashdagi noaniqlik va boshqalar).

Diskret sxemaning aniqligi haqida gapirganda, ular odatda hosilalarni 27 ga yaqinlashtirishdagi xatoni anglatadi. Xususan, agar yaqinlashish xatosi hisoblash to'r qadamining ikkinchi kuchi bilan taqqoslanadigan bo'lsa, u holda diskret sxema ikkinchi darajali aniqlikka ega deyiladi. Ushbu masala § 3da batafsilroq muhokama qilingan.

5.2 Muvofiqlik

Diskret sxema deyiladi kelishilgan asl differensial tenglama bilan, agar hisoblash tarmog'i aniqlanganda, yaqinlashish xatosi (3-bandga qarang) nolga moyil bo'lsa,

Muvofiqlikka erishish uchun qo'shimcha shartlar bajarilishi kerak bo'lgan ma'lum hisoblash sxemalari mavjud [Anderson va K]. Hisoblash sxemalarining izchilligini tekshirish dasturiy ta'minotni ishlab chiquvchilarning (va foydalanuvchilarning emas) vazifasi bo'lganligi sababli, bu masala bu erda batafsil ko'rib chiqilmaydi.


Yopish