Ribinio tūrio metodas ant netaisyklingo tinklelio. Ribinio tūrio metodas. Rekomenduojamas disertacijų sąrašas
Prieš kurį laiką ieškojau operacijų ir procesų, vykstančių OpenFOAM skaitmeninio modeliavimo bibliotekoje, aprašymo. Radau daug abstrakčių baigtinio tūrio metodo veikimo aprašymų, klasikinių skirtumų schemų, įvairių fizikinių lygčių. Norėjau sužinoti išsamiau - iš kur atsirado šios reikšmės tokiame ir tokiame išvesties faile esant tokiai ir tokiai iteracijai, kokios išraiškos yra už tam tikrų parametrų fvSchemes, fvSolution nustatymų failuose?
Tiems, kurie taip pat domisi - šis straipsnis. Tie, kurie gerai išmano OpenFOAM ar jame įdiegtus metodus – apie rastas klaidas ir netikslumus rašykite asmenine žinute.
Habré svetainėje jau buvo keli straipsniai apie OpenFOAM:
Todėl nesigilinsiu į tai, kad tai „atvira (GPL) skaitmeninio modeliavimo platforma, skirta modeliavimui, susijusiam su dalinių diferencialinių lygčių sprendimu baigtinio tūrio metodu, ir yra plačiai naudojama sprendžiant kontinuumo mechanikos problemas“.
Šiandien aš naudosiu paprastą pavyzdį, kad apibūdinčiau operacijas, kurios vyksta atliekant skaičiavimus OpenFOAM.
Taigi, atsižvelgiant į geometriją - kubas, kurio kraštinė yra 1 metras:
Mes susiduriame su užduotimi modeliuoti tam tikro skaliarinio lauko (temperatūros, medžiagos kiekio) srauto sklidimą, kurį pateikia tokia transporto lygtis (1) kūno tūrio viduje.
(1) 
,
Kai skaliarinis dydis, pavyzdžiui, išreiškia temperatūrą [K] arba tam tikros medžiagos koncentraciją ir išreiškia medžiagos perdavimą, masės srautą [kg/s].
Ši lygtis, pavyzdžiui, naudojama šilumos sklidimui modeliuoti 
,
kur k yra šilumos laidumas ir temperatūra [K].
Divergencijos operatorius iš tikrųjų yra
operatorius .
Leiskite jums priminti, kad yra nabla operatorius (Hamilton operatorius), kuris parašytas taip: 
,
Kur i, j, k yra vienetiniai vektoriai.
Jei nabla operatorių skaliariai padauginsime iš vektoriaus dydžio, gausime šio vektoriaus divergenciją:
„Fizikos požiūriu vektorinio lauko divergencija yra rodiklis, kiek tam tikras erdvės taškas yra šio lauko šaltinis arba grimzlė“.
Jei nabla operatorių padauginsite iš skaliro, gausite to skaliro gradientą:
Gradientas rodo skaliarinio dydžio padidėjimą arba sumažėjimą tam tikra kryptimi.
Problemos ribinės sąlygos yra tokios: yra įvesties paviršius, išėjimo paviršius, o likę paviršiai yra lygios sienos.
Kubo tūrio padalijimas į baigtinius tūrius
Mūsų tinklelis bus labai paprastas – kubą padaliname į 5 lygias ląsteles išilgai Z ašies.
Daug formulių
Baigtinio tūrio metodas numato, kad (1) integralioje formoje (2) bus patenkintas kiekvienam baigtiniam tūriui.
(2)
,
Kur yra galutinio tūrio geometrinis centras.
Galutinio tūrio centras
Supaprastinkime ir paverskime pirmąjį išraiškos (2) narį taip:
(2.1) (HJ-3.12)*
Kaip matote, padarėme prielaidą, kad skaliarinis dydis baigtinio tūrio viduje kinta tiesiškai, o kiekio vertę tam tikru tašku baigtinio tūrio viduje galima apskaičiuoti taip:
Norėdami supaprastinti antrąjį išraiškos narį (2), naudojame apibendrintą Gauso-Ostrogradskio teoremą: vektoriaus lauko divergencijos per tūrį integralas yra lygus vektoriaus srautui per paviršių, ribojantį duotą tūrį. Žmonių kalba, „visų srautų į/iš baigtinio tūrio suma yra lygi srautų per šio baigtinio tūrio paviršius sumai“:
(2.3)
,
Kur yra uždaras paviršius, ribojantis tūrį,
- vektorius, nukreiptas išilgai normalios nuo tūrio.
Vektorius S
Atsižvelgiant į tai, kad baigtinį tūrį riboja plokščių paviršių aibė, išraiška (2.3) gali būti transformuota į viso paviršiaus integralų sumą:
(2.4) (HJ-3.13)
,
Kur išreiškia kintamojo vertę veido centre,
- ploto vektorius, išeinantis iš veido centro, nukreiptas nuo ląstelės (lokaliai), nuo mažesnės indekso ląstelės į aukštesnio indekso ląstelę (globalus).
Šiek tiek daugiau apie vektorių S
Kad nebūtų du kartus išsaugoti tie patys vektoriaus parametrai, nes Akivaizdu, kad dviejų gretimų ląstelių normalus vektorius iki krašto tarp ląstelių, nukreiptas nuo ląstelės centro, skirsis tik krypties ženklu. Todėl tarp krašto ir ląstelės buvo sukurti savininko ir kaimyno santykiai. Jei ploto vektorius (globali, teigiama kryptis iš ląstelės su mažesniu indeksu į ląstelę su didesniu indeksu) rodo IŠ ląstelės centro, toks ryšys tarp ląstelės ir vektoriaus, o tiksliau tarp ląstelės ir veidas, žymimas savininkas). Jei šis vektorius rodo atitinkamos ląstelės viduje, tada kaimynas. Kryptis turi įtakos vertės ženklui (+ savininkui ir - kaimynui) ir tai svarbu sumuojant, žr. toliau.
Apie skirtumų schemas
Vertė veido centre apskaičiuojama pagal gretimų ląstelių centrų reikšmes - šis išraiškos metodas vadinamas skirtumo schema. OpenFOAM skirtumų schemos tipas nurodytas faile /system/fvSchemes:DivSchemes (numatytasis nėra; div(phi,psi) Gauso tiesinis;)
Gausas- reiškia, kad pasirinkta centrinio skirtumo schema;
linijinis- reiškia, kad interpoliacija iš ląstelių centrų į veidų centrus vyks tiesiškai.
Tarkime, kad mūsų skaliarinis dydis kinta tiesiškai baigtinio tūrio viduje nuo centro iki kraštų. Tada apytikslė vertė veido centre bus apskaičiuojama pagal formulę:
Kur yra svoriai ir apskaičiuojami kaip
Kur yra ląstelių tūriai.
Iškreiptų langelių atvejais yra sudėtingesnės formulės apytiksliesiems svoriams apskaičiuoti.
Taigi, phi_f reikšmės langelio kraštų centruose apskaičiuojamos pagal vertes ląstelių centruose. Gradiento reikšmės grad(phi) apskaičiuojamos remiantis phi_f reikšmėmis.
Ir visą šį algoritmą galima pavaizduoti šio pseudokodo forma.
1. Deklaruojame baigtinių tūrių gradientų masyvą, inicijuojame jį nuliais 2. Einame per visus vidinius paviršius (kurie nėra ribos) > Apskaičiuojame flux_f = phi_f*S_f. Apskaičiuokite phi_f reikšmes pagal phi reikšmes langelio centais > Pridėkite flux_f prie savininko elemento gradiento ir -flux_f prie kaimyninio elemento gradiento 3. Pakartokite per visus kraštinius > Apskaičiuoti flux_f = phi_f*S_f > Pridėkite flux_f prie savininko elemento gradiento (kaimynas -ribiniai paviršiai neturi elementų) 4. Pereikime visus elementus > Gautą gradiento sumą padalinkite iš elemento tūrio
Laiko atranka
Atsižvelgiant į (2.1) ir (2.4), (2) išraiška yra tokia:
(3)
Pagal baigtinio tūrio metodą atliekama laiko diskretizacija ir išraiška (3) rašoma taip:
(4)
Integruokime (4):
(4.1)
Padalinkime kairę ir dešinę puses į:
(5)
Duomenys imties matricai
Dabar galime gauti tiesinių lygčių sistemą kiekvienam baigtiniam tūriui.Žemiau yra tinklelio mazgų, kuriuos naudosime, numeracija. 
Mazgo koordinatės saugomos /constant/polyMesh/points
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))
Ląstelių centrų mazgų numeracija (50, 51 – kraštinių paviršių centrai): 
Veido centro mazgų numeracija: 
Elementų tūriai:
Interpoliacijos koeficientai, reikalingi langelių paviršių reikšmėms apskaičiuoti. Indeksas „e“ reiškia „dešinį langelio kraštą“. Dešinysis rodinio atžvilgiu, kaip parodyta paveikslėlyje „Ląstelių mazgų centrų numeracija“:
Atrankos matricos formavimas
Jei P = 0.
Išraiška (5), apibūdinanti kiekio elgesį
Bus transformuota į linijinių algebrinių lygčių sistemą, kurios kiekviena iš šių formų:
Arba pagal taškų indeksus ant veidų
Ir visi srautai į / iš ląstelės gali būti išreikšti suma
Kur, pavyzdžiui, srauto tiesinimo koeficientas langelio E centre,
- srauto tiesinimo koeficientas veido centre,
- netiesinė dalis (pavyzdžiui, konstanta).
Pagal veidų numeraciją išraiška bus tokia: 
Atsižvelgiant į elemento P_0 ribines sąlygas, tiesinė algebrinė lygtis gali būti pavaizduota kaip
...pakeisti anksčiau gautus koeficientus...
Srautas iš įleidimo angos"a nukreipiamas į ląstelę, todėl turi neigiamą ženklą.
Kadangi mūsų kontrolinė išraiška, be difuzijos termino, turime ir laiko terminą, tačiau galutinė lygtis atrodo taip
Jei P = 1.
Jei P = 4.
Tiesinių algebrinių lygčių sistema (SLAE) gali būti pavaizduota matricos forma kaip
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 = matmenys; internalField nevienodas sąrašas
Remiantis tuo, gaunamos vektoriaus reikšmės
Tada vektorius pakeičiamas SLAE ir įvyksta nauja vektoriaus skaičiavimo iteracija.
Ir taip toliau, kol neatitikimas pasiekia reikiamas ribas.
Nuorodos
* Kai kurios šiame straipsnyje pateiktos lygtys yra paimtos iš Jasak Hrvoje disertacijos (HJ yra lygties numeris) ir, jei kas nori daugiau apie jas paskaityti (Anksčiau buvo minimas subdomeno metodas, kuris buvo daugelio skaitmeninių metodų atskaitos taškas. Vienas iš tokių metodų yra baigtinio tūrio metodas. Tas pats metodas yra kitos plačiai paplitusios klasės – integralinių metodų – atstovas. Iš klasikinės subdomeno metodo žymėjimo formos imamas skaičiavimo srities padalijimas į subdomenus ir likučio integravimas per subdomeną. Skirtumas yra tai, kad nėra aiškaus aproksimacinės (bandymo) funkcijos įrašo. Tačiau, kaip ir anksčiau, mes stengiamės „tiksliai“ išspręsti lygtį kiekviename padomenyje. Todėl pradinė lygtis yra integruota per padomenį. Integraliniams metodams būdinga tai, kad pirmiausia imamas diferencialinės lygties integralas ir gaunama integralinė lygties užrašymo forma. Tada šios formos lygtis pritaikoma atskiroms tinklelio ląstelėms. Šiuo atveju ląstelės ir posritys yra tas pats.
Tiesą sakant, integralinė lygčių rašymo forma (fizikos požiūriu) yra dar platesnė nei diferencinė. Faktas yra tas, kad esant funkcijų netolygumams, diferencialinės lygtys netaikomos, o jų vientisieji analogai ir toliau veikia, veikia ir veikia.... Deja, juos įgyvendinus skaitmeniniu būdu, šis pranašumas kartais prarandamas.
Paprastai lygčių integralai turi paprastą ir suprantamą fizinę reikšmę. Pavyzdžiui, apsvarstykite tęstinumo lygtį. Parašyta pradinė diferencialinė lygtis
Integruokime jį per tūrį V, kurio paviršius S, ir laikui bėgant intervale nuo t 0 iki t 1. Integruodami išvestines naudojame Stokso formulę (ypatingi jos atvejai vadinami Green ir Ostrogradsky-Gauss formulėmis). Kaip rezultatas, mes gauname
Šiame žymėjime skirtumas tarp pirmųjų dviejų integralų reiškia masės pokytį tam tikrame tūryje per nagrinėjamą laiko intervalą. O dvigubas integralas parodo masę, per tą patį laikotarpį tekančią į tam tikrą tūrį per jį ribojantį paviršių. Natūralu, kad kadangi kalbame apie skaitinius metodus, šie integralai apskaičiuojami apytiksliai. Ir čia prasideda aproksimacijos klausimai, panašūs į tuos, kurie nagrinėjami baigtinio skirtumo metodu.
Panagrinėkime vieną iš paprasčiausių atvejų – dvimatę stačiakampę vienodą tinklelį. Taikant baigtinio tūrio metodą, funkcijų reikšmės paprastai nustatomos ne tinklelio mazguose, o langelių centruose. Atitinkamai taip pat indeksuojamos ne tinklelio linijos kiekviena kryptimi, o langelių sluoksniai (žr. pav.).
| j-1 |
| j |
| j+1 |
| k-1 |
| k |
| k+1 |
| A |
| B |
| C |
| D |
Šiuo atveju lygties integralioji forma bus parašyta taip
Kaip matote, šiuo atveju gavome įprastą lygtį, kurią taip pat galėjome parašyti baigtinio skirtumo metodu. Tai reiškia, kad jai gali būti taikomi tie patys stabilumo tyrimo metodai. (Greitas klausimas: ar ši schema stabili?)
Bet jei gautume tą patį, ar buvo verta statyti visą šį sodą? Paprasčiausiais atvejais naudos tikrai negauname. Tačiau sudėtingesnėse situacijose išryškėja nauda. Pirma, kaip minėta pirmiau, tokie metodai (net ir tokiu paprastu įgyvendinimu) daug geriau apibūdina nenutrūkstamumą ir sritis su dideliais nuolydžiais. Tuo pačiu metu garantuojamas masės, impulso ir energijos išsaugojimo dėsnių įvykdymas, nes jie yra stebimi kiekvienoje ląstelėje. Antra, šie metodai gali atlaikyti įvairius piktnaudžiavimus tinkle. Netgi kreivinės, nelygios ir netaisyklingos tinkleliai neišmuša šių metodų iš vėžių. Ši nauda ypač dažnai jaučiama, kai nurodomos ribinės sąlygos.
| j-1 |
| j |
| j+1 |
| k-1 |
| k |
| k+1 |
| A |
| B |
| C |
| D |
| E |
Pavyzdžiui, paveiksle pavaizduotu atveju integrali lygties forma turės formą
tai yra tiesiog ten, kur mes perėmėme integralą per visos ląstelės plotą, dabar perimame jį per „apkarpytą“ sritį, kur integralą perėmėme per visą kraštą, dabar perimame jį per likusią jo dalį . Buvo pridėtas integralas per ribą. Bet tai nesunku rasti pagal ribines sąlygas. Visų pirma, jei per sieną nėra tiekiamas masės srautas (o taip pat jokia masė nenunešama nuo paviršiaus ir (arba) neatsižvelgiame į jonų, prarandančių krūvį sienoje, masinį srautą), toks integralas yra tiesiog lygus nuliui. Panašioje energijos lygties formoje, kaip taisyklė, reikia atsižvelgti į srautą per sieną. Bet taip pat nesunku rasti iš ribinių sąlygų (jei jos nustatytos teisingai).
Norėdami tai sustiprinti, apibūdinkime, kaip atrodys baigtinio tūrio metodo taikymas vienai iš impulso išsaugojimo lygčių. Paimkime plokščią stacionarų atskirai įkrautų jonų korpusą. Neatsižvelgiame į klampumą ir elastingus susidūrimus. Gauname lygtį
Dėl stačiakampio tinklelio (žr. paveikslėlį aukščiau) gauname
Paprasčiausią tokios lygties aproksimaciją galima parašyti taip:
po sumažinimų gauname formulę
Naudojimas baigtinio (valdymo) tūrio metodas Parodykime dvimatės stacionarios šilumos lygties pavyzdžiu:
Ryžiai. 13. Skaičiavimo tinklelis, naudojamas (31) lygčiai išspręsti
baigtinio tūrio metodas
Naudodami vidutinės reikšmės teoremą galime parašyti
|
|
,kur Δx, Δу yra langelio paviršių ilgiai, x W yra A langelio kairiosios („vakarinės“) kraštinės abscisė, x E yra dešinės („rytinės“) kraštinės abscisė, y N yra ordinatės viršutinės („šiaurinės“) ribos, y S yra apatinės („pietinės“) ribos ordinatė, S * – vidutinis ląstelės šilumos išsiskyrimo greitis. Išvestinių priemonių indeksas (*), esantis kairėje (32) pusėje, rodo, kad jie turėtų būti laikomi vidutinėmis vertėmis, nustatytomis taip, kad teisingai atspindėtų šilumos srautus kiekvienoje iš ribos. Atsižvelgiant į šią aplinkybę, atskirą (32) analogą galima gauti be vargo [Patankar].
Taigi (32) lygtis apibūdina šilumos balansą (energijos tvermės dėsnį) ląstelėje A. Jei šilumos srautai tarp elementų yra teisingai aprašyti, sistema, sudaryta iš (32) formos lygčių, pritaikytų kiekvienam kontroliniam tūriui apibūdinkite šilumos balansą visoje skaičiavimo srityje.
Pastraipos pabaigoje reikia pažymėti, kad tam tikrais atvejais aukščiau aprašytais metodais gautos skaičiavimo formulės gali sutapti, o reikšmingiausi skirtumai atsiranda naudojant kreivinius nestačiakampius skaičiavimo tinklelius.
5. Diskrečiųjų grandinių savybės
5.1 Tikslumas
Tikslumas apibūdina skaitinės schemos priimtinumą jos praktiniam naudojimui. Atrodo, kad diskrečios grandinės tikslumo įvertinimas yra labai sudėtinga užduotis, nes paaiškėja, kad beveik neįmanoma atskirti klaidų, atsirandančių dėl grandinės savybių, nuo klaidų, atsirandančių dėl kitų veiksnių (pvz. apvalinimo paklaidos, netikslumas nurodant ribines ir pradines sąlygas ir pan.).
Kalbėdami apie diskrečios schemos tikslumą, dažniausiai jie turi omenyje išvestinių aproksimacijos paklaidą 27 . Visų pirma, jei aproksimavimo paklaida yra palyginama su antrąja skaičiavimo tinklelio žingsnio galia, tada sakoma, kad diskrečioji schema turi antros eilės tikslumą. Išsamiau šis klausimas buvo aptartas § 3.
5.2 Nuoseklumas
Diskretinė grandinė vadinama susitarta naudojant pradinę diferencialinę lygtį, jei patikslinus skaičiavimo tinklelį aproksimacijos paklaida (žr. § 3) yra lygi nuliui,
Yra žinomos skaičiavimo schemos, kuriose turi būti įvykdytos papildomos sąlygos, kad būtų pasiektas nuoseklumas [Anderson ir K]. Kadangi skaičiavimo schemų nuoseklumo tikrinimas yra programinės įrangos kūrėjų (o ne vartotojų) užduotis, šis klausimas čia plačiau nebus nagrinėjamas.
