Dílčí projekt: Systém projektování textilních struktur 3. Vývojové etapy VÝZKUMNÉ CENTRUM TEXTIL Textilní materiály a konstrukce textilních výrobků Dílčí projekt: Systém projektování textilních struktur ZÁVĚREČNÁ ZPRÁVA Milda Pociute, Jiří Chaloupek, Eva Košťáková, Larisa Očeretna, David Lukáš 3.6 RAYLEIGHOVA NESTABILITA KAPALINOVÝCH FILMŮ POKRÝVAJÍCÍCH VLÁKNO 3.6.1. ÚVOD Vzájemné působení kapalin nebo polymerních tavenin a vlákenných systémů je důležité v mnoha technologických procesech jako jsou: barvení a zušlechťování textilií, lakování drátů, výroba pojených textilií, atd.. Neméně významná je tato interakce i v řadě případů používání hotových výrobků. Jako příklad jmenujme: dětské pleny, dámskou hygienu, vlasovou kosmetiku, kapalinovou filtraci. Je s podivem, že přes velké množství výrobních i uživatelských situací, ve kterých se smáčení vláken uplatňuje, je teoretický popis smáčení vlákenných útvarů stále v ranném stádiu vývoje [1]. Hlavní zvláštností odlišující jevy smáčení rovinného povrchu a válce (který je geometricky nejjednodušším modelem vlákna) jsou následující . Když pokryjeme vlákno nebo drát tenkým, ale stále ještě makroskopickým, kapalinovým filmem, začne se tento samovolně rozpadat na malé kapkovité útvary, které jsou na vlákně rozmístěny více méně v pravidelných rozestupech. Vytváření kápek z původně rovnoměrného kapalinového filmu nastane i v případě, že úhel smáčení mezi kapalinou a vláknem je nulový, a tedy takový, který by v případě smáčení rovinného povrchu vyvolal dokonalé smáčení [2]. Výše popsané, na první pohled neobvyklé, smáčení vláken je ve skutečnosti projevem nestability kapalinových sloupců popsané prvně Plateauem [3] a Rayleighem [4]. Kapalinový sloupec o poloměru r, na který působí porucha o vlnové délce λ, je nestabilní a v důsledku povrchových sil se rozpadne na kapičky v případě, že vzdálenost středů kapek λ je větší než 2πr [2]. Problém rozpadu kapalinového sloupce byl studován mnoha autory jak z teoretického tak experimentálního hlediska.Uvedeme zde aspoň dvě takové práce: Tomotikovu [5] a Meisterovu [7]. My budeme sledovat mírně pozměněný problém, kde uvnitř kapalinového sloupce bude jako pevné jádro umístěno vlákno. Tato problematika byla v sedmdesátých letech minulého století studována analytickými prostředky Roem [2]. V této práci bude ke studiu stejné problematiky využita počítačová simulace založena na metodě Monte-Carlo a Auto-modelu. Cílem této práce je rozvinout novou metodu pro popis nestabilit a morfologických přechodů kapalinových těles v interakci se strukturně mnohotvárnými nebo heterogenními substráty. Tato metoda je založena na uplatnění postupů Monte Carlo ve statistické fyzice s využitím auto-modelů uplatňovaných v teorii rozpoznávání textur. Ukážeme, že Hammersley Cliffordova věta je přímým důsledkem Besagova tvrzení. K tomu rozložíme příslušná Markovova pole na carové systémy o stejném počtu interagujících uzlů a pomocí jednoduché grafické metody ukážeme různou roli „ heterogenních“ a „ homogenních“ klanů. Upřesníme tak třídu auto-modelů, která je bohatší než běžně užívaný Isingův model ve fyzice pevné fáze. Cílem příspěvku je také ukázat, že zvolený auto-model dokáže správně
Dílčí projekt: Systém projektování textilních struktur 3. Vývojové etapy reflektovat úhrnou bilanci přebytku energie ve fázové vrstvě rozhraní mezi kapalinou a plynem při rozpadu kapalinového sloupce na jednotlivé kapky tím, že porovnáme Rayleighovu vlnovou délku s délkou získanou simulačním experimentem a zjistíme podivuhodnou shodu. Naše práce je rozdělena do následujících oddílů: Nejprve uvedeme úvahy o rotačně symetrických plochách o konstantní hodnotě střední křivosti, které jsou stabilní. Následovat budou hrubé úvahy o nestabilitě kapalinových sloupců založené na porovnání energií a entropie kapalinového sloupce a kapek, které vzniknou jeho rozpadem. V dalším oddíle bude podrobněji popsána Rayleighem vypracovaná metoda pro předpověď charakteristické délky (vlnové délky), která řídí rozpad kapalinového sloupce na kapky. Dále se budeme zabývat simulační metodou pro modelování nestability kapalinových sloupců a její výsledky budou srovnány s klasickými výstupy Rayleighovy a Royovy teorie. 3.6.2. ROTAČNĚ SYMETRICKĚ PLOCHY O KONSTATNÍ HODNOTĚ KŘIVOSTI Vyloučíme-li z našich úvah vliv vnějších polí (například gravitačního pole), které by mohly ovlivnit tvar kapalinového filmu na vlákně, pak je přirozené očekávat, že rovnovážné kapalinové těleso přisedlé k vláknu bude rotačně symetrické. Podmínkou jeho rovnováhy pak bude stejná hodnota kapilárního tlaku v celém tělese. Hodnota kapilárního tlaku Pc v závislosti na hodnotách hlavních křivostí K1 a K2 (nebo jím příslušejících poloměrů R1 a R2 oskulačních kružnic) byla poprvé odvozena v pracích Younga [8] a Laplace [9]. ⎛ 1 1 ⎞ ⎟⎟ . (1) Pc = γ ⎜⎜ + ⎝ R1 R2 ⎠ Okamžitým důsledkem výše uvedeného vztahu pro mechanickou rovnováhu K + K2 1 kapalinového tělesa je konstantní hodnota střední křivosti K = 1 , kde K 1 = a 2 R1 1 K2 = , platící pro celý jeho povrch. R2 Rotačně symetrické plochy o konstantní hodnotě střední křivosti jsou označovány jako H - rotační plochy (H- surfaces of revolution) [10]. Plateau ukázal, že existuje právě šest různých druhů H-rotačních ploch: rovina a katenoid (střední křivost je rovna nule), koule válec unduloid a nodoid (střední křivost je nenulová). V roce 1841 ukázal francouzský matematik Delanuy možnost vytvářet všechny H-rotační plochy velmi jednoduchým způsobem, který je pospán následovně: vybereme si přímku, která bude dále sloužit jako osa rotace. Dále si zvolme některou z kuželoseček - to znamená elipsu (která může být degenerována na úsečku), kruh, parabolu nebo hyperbolu, kterou necháme po dříve zvolené ose rotace odvalovat. Ohniska výše zmíněných kuželoseček odvalovaných po přímce vytváří křivku, která si říká "roulade". Když tuto křivku roztočíme vzhledem k dříve zvolené ose rotace, vytvoří tato křivka H – rotační plochu. Pět křivek typu "roulade" ukazuje Obr.1.
Dílčí projekt: Systém projektování textilních struktur 3. Vývojové etapy Válec Kruh
Unduloid Elipsa
Koule Úsečka
Katenoid Parabola
Nodoid Hyperbola
Obr.1.: Obrázek znázorňuje křivky "roulade" a jim příslušející názvy H-rotačních ploch.
Například ohnisko kruhu vytvoří "roulade" ve tvaru rovnoběžky s osou rotace. Rotací této rovnoběžky kolem osy rotace získáme válec. Odvalující se ohnisko elipsy vytvoří "roulade" ve tvaru horního profilu řady kapek přisedlých na vlákně. Roztočením této křivky vzhledem k ose rotace obdržíme unduloid. Další H-rotační plocha blízká problematice této práce se obdrží krácením úsečky (degenoravané elipsy jejíž ohniska jsou v krajních bodech úsečky) po přímce. Následnou rotací získáme plochy ve tvaru povrchů koulí. Ve stavu beztíže by kapalina vzlínající z nádrže na vlákno měla podobu katenoidu [11]. Katenoid získáme rotací "roulade", která je stopou odvalujícího se ohniska paraboly. Konečně postupným střídáním stop ohnisek odvalujících se dvou větví hyperboly získáme křivku jejíž rotací obdržíme nodoid. Z uvedeného je patrné, že změny tvarů (nestabilita kapalinového filmu na vlákně (se může odehrávat pouze v rámci uvedených šesti H-rotačních ploch, pokud považujeme za podmínku že výsledné kapalinové těleso je kvazistacionární (jen pomalu měnící svůj tvar v čase). 3.6.3 ENERGETICKÁ BILANCE RAYLEIGHOVY NESTABILITY Hrubé vysvětlení jevu, který je nazýván Rayleighova nestabilita kapalinových sloupců, může být provedeno sledováním hodnot energií počátečního tělesa - tedy válce – a výsledkem nestability, kterou je řada stejně velkých a rovnoměrně od sebe vzdálených koulí za podmínky zachování objemu kapaliny původního a výsledného systému. Pro objemy paltí následující rovnost 4 πrc2 λ = πrs3 , (2) 3
kde λ je úsek válcového kapalinového tělesa (vlnová délka) o poloměru rc, který se přemění právě na jednu kulovou kapku o poloměru rs. Povrchová a tlaková energie (poslední z nich se vyjádří jako součin kapilárního tlaku ρ 4 a objemu V) se vyjádří pomocí výrazů: 2πrc λγ , 4πrs2γ a γrc−1πrc2 λ , 2γrs−1 πrs3 . Porovnáním 3 součtu povrchové a tlakové energie válce a koule dostaneme
Dílčí projekt: Systém projektování textilních struktur 3. Vývojové etapy 8 2πrc λγ + γπrc λ ≥ 4πrs2γ + γπrs2 . 3
(3)
Z rovnice (2) vyjádříme: 2
⎛3 ⎞3 r ≤ ⎜ rc2 λ ⎟ . (4) ⎝4 ⎠ a dosadíme do upravené rovnice (3), která má tvar: 20 rs2 λ≥ . (5) 9 rc Po dosazení a další jednoduché úpravě obdržíme: 5 3.4 λ ≥ 4 rc ≅ 1,96πrc (7) 3 tento výsledek je kvalitativně srovnatelný s Rayleighovou hodnotou pro vlnovou délku λ. Jak uvidíme v dalším odstavci je správná hodnota pro vlnovou délku λ rovna 2,86 πrc, které námi odvozená nerovnost vyhovuje. Přesto musíme podtrhnout tu skutečnost, že správný výsledek pro vlnovou délku z Rayleighovy nestability rozpadu kapalinového sloupce nelze získat pouze použitím zákona zachování volné energie. Přeměna tvaru kapalinového tělesa, která je doprovázena změnou velikost jeho povrchu způsobuje i změny celkové hodnoty přebytku entropie ve fázovém rozhraní mezi kapalinou a plynem, který je na tento povrch tohoto rozhraní vázán [12]. 2 s
3.6.4. RAYLEIGHOVA NESTABILITA Rayleighova nestabilita kapalinových sloupců vzniká jako důsledek jemných poruch (kapilárního vlnění [13]), o kterých předpokládáme, že mají harmonický průběh s amplitudou, která se exponenciálně zvětšuje. Při šíření takovéto poruchy kapalinového sloupce dochází k přelévání povrchové energie na energii kinetickou, která je spojena s prouděním kapaliny při změně jejího tvaru, to je při přeměně válcového tělesa na řadu oddělených kapek. Při této přeměně tvaru se jednotlivé poruchy nevyvíjejí stejně rychle. Jejich vývoj je totiž závislý na vlnové délce šířící se poruchy. Porucha s nejrychlejším vývojem je potom ta, která určí výslednou vlnovou délku rozpadu kapalinového sloupce kapky. Zdůrazněme ještě ten fakt, že vítěznou poruchu stačí sledovat jen v jejích počátečních stádiích vývoje. Podrobnější popis Rayleighovy nestability je následující: Porucha se šíří tak, že poloměr r kapalinového tělesa se podél jeho osy rotace změní podle vztahu r = r0 + a. exp[qt ] cos(kz ) , (8) kde ro je poloměr neporušeného sloupce, a značí počáteční amplitudu poruchy, q je parametr určující exponenciální růst amplitudy s časem , k je vlnový vektor (k=2πλ-1). Pro jednoduchost budeme dále značit amplitudu poruchy v daném čase α (α = a. exp (q.t)). Protože celý proces sleduje v ranném stádiu vývoje, kdy je amplituda malá α << 1, budeme v dalším zanedbávat členy řádu α2. Pro změnu povrchové energie V kapalinové trysky s poruchou o amplitudě α pak získáme vztah V = +γ
πα 2
(1 − k 2 r02 ). (9) 2r0 Kinetickou energii T odhadneme z pohybu tenké kapalinové trubice o tlouštce dr a poloměru ro , který má na délku λ hmotnost m a pohybuje se kolmo na svoji osu rotace rychlostí v(z). Tuto rychlost zjistíme časovou derivací hodnoty r(t).
Dílčí projekt: Systém projektování textilních struktur 3. Vývojové etapy dr dα =´α 2 cos 2 (kz )dz , ´α = = q 2α 2 . (10) dt dt Pro kinetickou energii pak platí 1 1 T = πρr0 dr q 2α 2 (11) 2 k Porovnáním úbytku potencionální a kinetické energie dostaneme vztah pro parametr rychlosti vývoje amplitu g v závislosti na součinu vlnového vektoru z k a poloměru ro . v=
q=
γ ρr dr 3 0
. kr0 . 1 − k 2 r02 .
(12)
Hledání extrému q (v tomto případě maxima) z podmínky nulovosti derivace nás přivede k rovnici: d ( x(1 − x 2 ) = 0) , dx v niž jsme převzali kro za x. Výše uvedená rovnice má řešení x = 1 3 ≅ 0,577. Odtud získáme pro vlnovou délku λ hodnotu 3,46 πro. Oproti Rayleighovu výsledku λ = 2,88 πro je odchylka našeho výpočtu přibližně +20%. Příčina této chyby tkví v našem hrubém odhadu kinetické energie. Správný výsledek bychom obdrželi použitím rychlostního pole proudění kapaliny, které by bylo získáno řešením Navier-Stokesovy rovnice při daných okrajových podmínkách [14]. 3.6.5 AUTO-MODELY Pro počítačovou simulaci Rayleighovy nestability jsme použily mřížkový model, ve kterém je prostor i čas diskretizován. Simulační krabice (simulační box) obsahuje N uzlů uspořádaných do čtvercové trojrozměrné mříže o Nx,Ny, a Nz uzlech rozložených podél hran simulační krabice po řadě rovnoběžně s osami x, y a z. Množinu všech uzlů mříže budeme označovat S a jednotlivé uzly budou značeny s1, s2, …, sN. V každém uzlu s∈S je přiřazena hodnota mřížové proměnné Xs. Různým hodnotám mřížové proměnné v daném uzlu budeme později přisuzovat fyzikální význam: uzel může být vyplněn plynem, kapalinou, nebo pevnou látkou. Každé z těchto možností je přiřazen číselný kód 0, 1, nebo 2. Tyto tři možné hodnoty xs mřížové proměnné Xs tvoří tak zvaný lokální konfigurační prostor Λs={0, 1, 2}. G Konfiguraci mříže, to je množinu hodnot všech proměnných Xs, budeme označovat x , kde G x =(x1, x2,…,xn). Celkový konfigurační prostor systému bude dále značen Ω=ΛN. G Jednotlivé konfigurace systému x se mohou v souladu s principy statistické fyziky G vyskytovat s různými pravděpodobnostmi Π( x ). Právě popsaný systém uzlů s mřížovými proměnnými nabývajícími náhodných hodnot představuje tak zvané náhodné pole popsané například Gemanem v [14] nebo [15]. V polovině sedmdesátých let ukázal Besag [16], že pravděpodobnosti náhodných polí G Π( x ) jsou jednoznačně určeny svými lokálně podmíněnými pravděpodobnostmi Pi(xi/x1, G …,xi-1,xi+1,…,xN), které budeme zkráceně zapisovat také jako Pi(xi/ x (i)). Hodnoty lokálně podmíněné pravděpodobnosti pro uzel i udávají pravděpodobnost výskytu xi proměnné Xi při vybraných konstantních hodnotách mřížových proměnných ve všech ostatních uzlech. Jednoznačnost nalezená Besagem má podobu rovnice (14).
Dílčí projekt: Systém projektování textilních struktur 3. Vývojové etapy G Π ( x ) N Pi ( xi / x1 ,..., xi −1 , yi +1 ,... y N ) , G =∏ Π ( y ) i =1 Pi ( yi / x1 ,..., xi −1 , yi +1 ,... y N )
(14)
kde xs a ys jsou různé hodnoty mřížové proměnné v uzlu s. Důkaz tohoto tvrzení nalezne čtenář v [15, 16]. Speciálním případem náhodného pole je Markovovo náhodné pole. Důležitou vlastností Markovova náhodného pole je ta skutečnost, že pravděpodobnost se kterou mřížová proměnná Xs nabude hodnoty xs je podmíněna jen hodnotami xr náhodných proměnných Xr v sousedství uzlu s. Množinu sousedních uzlů pro uzel s budeme označovat Ns. Ze všech možných typů sousedství vybereme pro náš model sousedství symetrická; to je taková pro která platí s∈Nt⇔t∈Ns. Navíc nebudeme pokládat daný uzel za souseda sebe sama; s∉Ns. Takto vybraný typ sousedství a náhodného pole vyhovuje fyzikálním představám o mezimolekulárních potenciálech, které jsou omezeného dosahu. Dále budeme pokračovat definicí řádu sousedství. V této práci budeme uvažovat stejný typ symetrického sousedství Ns0 uzlu s a řádu o jako v [14 a 15].
{
2
}
N s0 = r ∈ S : 0 < s − r ≤ o ,
(15)
kde s − r je euklidovská vzdálenost uzlů s, r ∈ S. řád sousedství je pak kvadrát vzdálenosti mezi nejvzdálenějším sousedem r uzlu s a pozicí uzlu s. Dvojrozměrná sousedství 1., 2. a 4. řádu jsou zobrazena na obr.2.
s
s
o=1
o=2
s
o=3
Obr.2: Zleva doprava obrázek ukazuje postupně sousedství prvého, druhého a čtvrtého řádu. Uzel s je vždy vyznačen tmavě, aby se zdůraznilo, že daný uzel sám sobě není sousedem. Vybraný nevzdálenější soused je s uzlem s spojen úsečkou, kvadrát jejíž délky je roven řádu sousedství.
Dalším důležitým pojmem pro Markovova náhodná pole je klan [17] neboli klika podle [15]. Pro daný systém sousedství N={Ns⊂ S, s ∈ S} nazveme klanem takovou podmnožinu uzlů C z celkové jejich množiny S pro kterou platí, že každé dva uzly z C jsou si navzájem sousedy. Navíc definujeme, že jednotlivé uzly tvoří také klan. Uvažujeme-li dvojrozměrnou čtvercovou mříž a sousedství druhého řádu (o = 2), pak existuje celkem 10 různých typů klanů v takovém systému. Typ sousedství i typy klanů jsou uvedeny v obr.3.
o=2 Obr.3: Pro systém uzlů uspořádaných do dvourozměrné čtvercové mříže a pro symetrická sousedství druhého řádu (o=2) existuje 10 různých typů klanů.
Dílčí projekt: Systém projektování textilních struktur 3. Vývojové etapy Můžeme nyní shrnout předcházející uvedené poznatky v následující závěr: Lokální podmíněné pravděpodobnosti Pi(xi/x1,…,xi-1,xi+1,…,xN) pro obecná (Gemanova) náhodný pole jsou pro náhodná Markovova pole podmíněny jen stavy sousedů Pi(xi/xr∈Ni) a navíc typem klanů, který pro systém připustíme (klany tvořené jedním uzlem, dvojicí, trojicí, čtveřicí,….uzlů). Takováto podmínka znamená jistá omezení pro výběr pravděpodobností G Pi(xi/ x (i)) splňujících Besagovu podmínku (14). Toto omezení poprvé zkoumali v roce 1971 Hammersley a Clifford [18], tím že si položili následující otázku: Jaká je nejobecnější podoba podmíněných pravděpodobností Pi(xi/xr∈Ni) pro Markovo náhodné pole s vybraným typem symetrického sousedství? Protože naším předmětem zájmu budou auto-modely, což jsou modely připouštějící pouze klany jednoho a dvojic uzlů, omezíme náš výklad HammersleyCliffordovy věty na tuto třídu Markových polí. Hammerley-Cliffordovu větu pro zmíněný typ Markových polí zde představíme jako důsledek řešení Besagova vztahu (14). Nejprve omezme naši úvahu na symetrická sousedství nultého řádu, v kterém existují pouze klany jednoho uzlu. Besagovo tvrzení (14) tak nabude jednoduché podoby G N Pi ( x i ) Π i (x ) = G ∏ . (16) Π i ( y ) i =1 Pi ( y i ) G Tato funkcionální rovnice nabízí řešení Π ( x ) = P1 ( x1 ) ⋅ P2 ( x2 ) ⋅ ⋅ ⋅ ⋅ ⋅ PN ( x N ), které je užitečné pro další úvahy zapsat v podobě sumy N G Q ( x ) = ∑ g i (xi ), (17) i =1
G G kde Q( x ) = ln Π (x ) a gi(xi)=lnPi(xi) pro tuto úpravu se předpokládá, že Pi(xi)≠0 pro každý uzel i ∈ S. Obecným řešením Besagova vztahu (14) je tedy konstatování, že pravděpodobnosti G G Π (x ) výskytu dané konfigurace x zvoleného Markovova pole můžeme zapsat jako součin všech lokálně podmíněných pravděpodobností (v tomto případě – kdy uvažujeme pouze klany jednoho uzlu podmínka dokonce vymizela a pravděpodobnost lokální konfigurace je G nezávislá na stavu jakéhokoliv uzlu pole krom stavu uzlu sledovaného). Logaritmus Q( x ) pravděpodobnosti takové konfigurace pak lze zapsat jako sumu logaritmů lokálně podmíněných pravděpodobností. Zároveň považujeme nenulovost každé lokální pravděpodobnosti Pi(xi). Na funkce gi(xi) nejsou kladena žádná další omezení. Je však pro další výklad užitečné zapsat je v následujícím tvaru g i (xi ) = xi Gi ( xi ) . (18)
Funkce gi(xi) tak zapíšeme v podobě součinu hodnoty mřížové proměnné Xi v uzlu i a g (x ) funkce Gi(xi), pro kterou platí Gi (xi ) = i i . Řešení (17) pak nabývá podoby xi N G Q ( x ) = ∑ xi Gi ( xi ) i =1
(19)
Dílčí projekt: Systém projektování textilních struktur 3. Vývojové etapy Nyní přistoupíme k řešení Besagovy relace (14) pro Markovova náhodná pole se sousedstvím prvního řádu a klany tvořenými pouze dvojicemi uzlů. Besagův vztah (14) pro tento systém nabývá tvaru G P (x ) N Pi xi / xy(i ) , G =∏ (20) P ( y ) i =1 Pi yi / xy(i )
( (
) )
Kde jsme symbolem xy (i ) označily konfiguraci okolí uzlu i, která je dána stavovým vektorem G složeným z podmnožiny stavového vektoru x(i ) celého systému neobsahující prvek i G ( x(i ) = ( x1 ,..., xi −1 , yi +1 ,..., y N ) ), pro které uzly patří do okolí Ni prvku i. Rovnice (20) můžeme vyjádřit symbolicky pro systém devíti uzlů uspořádaných do dvojrozměrné čtvercové mříže. Toto symbolické vyjádření je zaznamenáno v obr.4.
Obr.4: Jednotlivá schémata v obrázku zobrazují konfigurace systému (x1,…,xi-1,xi,yi+1,…,yn), která jsou G G složena z částí konfigurací x a y . Úsečkami spojujícími uzel i s uzly jeho okolí Ni je zvýrazněna podmínka příslušné pravděpodobnosti. Například P4(x4/x5,y5,y7) značí pravděpodobnost výskytu stavu x4 uzlu 4 za podmínky, že jeho okolí N4 tvořené uzly 1, 5 a8 je v konfiguraci (x1,y5,y7). Značka x znázorňuje násobení mezi jednotlivými lokálně podmíněnými pravděpodobnostmi Pi a symbol ÷ je
(
dělením mezi součinem všech Pi x i / xy (i )
) a součinem všech P (y / xy ), tak jak je zapsáno v i
i
(i )
rovnici (20).
Pečlivým prohlédnutím obr.4 shledáme, že každému klanu v heterogenní G konfiguraci (tj. jeden uzel má stav příslušející konfiguraci x a druhý stav příslušející G konfiguraci y ) přísluší stejný klad s opačným umístěním ve zlomku pravé strany rovnice (16). Pro lepší orientaci jsme tyto klany (vybrali jsme uzly 5 a 6) v obr.4 zvýraznily. Totéž se děje i s klany typu
,
,
.
Své protějšky v čitateli a jmenovateli nemají pouze homogenní klany ,
. Všechny homogenní klany typu
a
,
a
jsou navíc umístěny v čitateli a klany
a jsou naopak všechny lokalizovány ve jmenovateli. Pohlížíme-li na Besagovo tvrzení (14) jako na funkcionální rovnici pro Markovova náhodná pole se symetrickými okolími prvního řádu a zároveň připouštíme pouze existenci klanů dvou prvků, pak je zcela přirozené navrhnout řešení takové relace ( tj. navrhnout tvar G pravděpodobnosti P(x ) ) v podobě součinu funkcí dvou proměnných G P( x ) = ∏ Pi , j (xi , x j ), kde uzly i a j tvoří klany. V součinu se každý klan vyskytne pouze typů
1≤i ≤ j ≤ N
Dílčí projekt: Systém projektování textilních struktur 3. Vývojové etapy jedenkrát. Měli bychom ještě zdůraznit tu skutečnost, že navrhnuté řešení je platné (má stejnou funkcionální podobu) pro libovolná očíslování uzlů Markovova náhodného pole. Nyní se pokusíme interpretovat význam funkcí Pi,j(xi, xj). za tímto účelem přejdeme od G systému uzlů k systému klanů. Každou konfiguraci x (x1 ,..., x N ) můžeme jednoznačně popsat i jako konfiguraci klanů, v našem případě systém ve kterém existují pouze klany dvou uzlů jako systém vazeb mezi dvojicemi uzlů. Vazby jako objekty vlastně představují také Markovo náhodné pole, a proto pro ně musí platit stejná pravidla jako pro uzly. Vzájemně se budou podmiňovat jenom ty klany, které obsahují společné uzly. Záměna „uzlových“ Markových polí za „klanové“ Markovy pole připomíná postupy užívané v teorii perkolace [18]. Očíslujeme-li všechny různé klany dvou prvků našeho systému, indexem α=1,…,θ, kde θ je celkový počet klanů dvou prvků našeho systému, pak pro ně musí platit také Besagovo tvrzení (14) G P (γ ) θ Pβ γ β / γσ ( β ) G =∏ , (21) Pδ β =1 Pβ δ β / γσ ( β )
()
( (
) )
G G kde γ a δ jsou dvě různé konfigurace klanů dvou prvků, β je index zvoleného klanu a γδ ( β ) G je konfigurace klanů z okolí klanu β. Je zřejmé, že konfiguraci klanů γ je jednoznačně G G G G přiřazena konfigurace uzlů x a musí platit P ( x ) = P (γ ) . Tak je naznačeno, že P ( x ) lze zapsat alternativně jako součin funkcí Pi,j(xi,xj) dvou proměnných, které mají význam pravděpodobnosti výskytu dané konfigurace (xi,xj) v klanu tvořeného dvojicí prvků (i,j) za podmínky jisté konfigurace sjednocených okolí uzlů i a j. Naznačili jsme, že každému Markovu poli s jistým systémem okolí a klanů různých typů přísluší Markovova pole jednotlivých typů klanů. Pravděpodobnost výskytu n-klanů ( kde n je počet uzlů v klanu) lze vyjádřit funkcemi n-proměnných. Těmito proměnnými jsou příslušné mřížové proměnné. Předpokládáme-li, že stavy polí klanů jsou na sobě v jistém G smyslu závislé, pak pravděpodobnost výskytu konfigurace x systému uzlů může být vyjádřena součinem právě takových klanových funkcí o stoupajícím počtu proměnných. G G Besagovo tvrzení (14) pak říká, že pravděpodobnost P ( x ) výskytu konfigurace x lze zapsat jako součiny součinů klanových funkcí ⎞ ⎞⎛ N G ⎛ N ⎞⎛ N P ( x ) ≅ ⎜⎜ ∏ Pi ( xi )⎟⎟⎜⎜ ∏ Pi , j (xi , x j )⎟⎟⎜⎜ ∏ Pi , j ,k (xi , x j , xk )⎟⎟ ⋅ ... (22) ⎠⎝ 1≤i ≤ j ≤ N ⎝ i =1 ⎠ ⎠⎝ 1≤i ≤ j ≤ k ≤ N Logaritmus tohoto vztahu (22) se zapíše jako N N G (23) Q( x ) = ∑ xi Gi ( xi ) + ∑ xi x j Gi , j (xi , x j ) + ∑ xi x j xk Gi , j ,k (xi , x j , xk ) + ... . i =1
1≤i ≤ j ≤ N
1≤i ≤ j ≤ k ≤ N
Definiční vztahy lnPi,j(xi,xj)=xixjGi,j(xi,xj), a podobně pro členy vyššího řádu jsme (22) přepsaly v (23), což je zápis obvyklý pro formulaci Hammersley-Cliffordovy věty. Hammersley Cliffordova věta může být vyslovena v následujícím znění: Pro každou skupinu indexů splňujících podmínku 1≤i≤j≤…≤s≤N, je funkce Gi,j,…,s nenulová tehdy, když uzly i,j,…,s tvoří klan. Krom tohoto omezení, mohou být G funkce vybírány libovolně. Libovůlí se zde myslí ta skutečnost, že každému výběru G funkcí přísluší jisté Markovo pole G s určitou pravděpodobností P ( x ) výskytu konfigurací. Markovu pole, ve kterém připouštíme jen klany jednoho a dvou uzlů, přísluší funkce G Q( x ) tvaru G Q( x ) = ∑ xi Gi (xi ) + ∑ xi x j Gi (xi , x j ). 1≤i ≤ N
1≤i ≤ j ≤ N
Dílčí projekt: Systém projektování textilních struktur 3. Vývojové etapy (24) Takové systémy pak budeme nazývat auto-modely. 3.6.6. POČÍTAČOVÁ SIMULACE RAYLEIGHOVY NESTABILITY POMOCÍ AUTO-MODELŮ Námi použitý auto-model pro studium Rayleighovy nestability se sestává ze simulační krabice o N x * N y * N z = 600 * 25 * 25 = 54 *10 4 uzlech. V každém uzlu i může mřížová proměnná X i nabývat jenom hodnoty 2, nebo alternativně hodnot 0,1. Hodnota 2 reprezentuje pevnou látku ( částici vlákna), která je během simulace nehybná. Pokud tedy je uzel obsazen hodnotou 2, pak je touto hodnotou obsazena trvale. Hodnota 1 představuje obsazení uzlu kapalinou a 0 obsazení uzlu vzduchem. G Funkce Q( x ) , kterou budeme dále nazývat hamiltonianem ( nebo-li celkovou energií systému) má podobu G 1 Q( x ) = (25) ∑ C ( xi , x j ) ,
τ 1≤i ≤ j ≤ N
kde C ( xi , x j ) jsou výměnné energie mezi sousedními uzly xi a x j . τ je termodynamická teplota. Hodnoty těchto výměnných energií jsou přehledně uvedeny v tabulce tab.1. Uvedený hamiltonian neuvažuje klany jednoho prvku, které by reprezentovaly externí pole ( například gravitační pole ). Protože z Rayleighova odvození vlnové délky rozpadu kapalinového sloupce plyne jeho nezávislost na materiálových parametrech kapaliny, vlnová délka je závislá podle vztahu λ = 2.88πr jen na poloměru sloupce, musí být výsledek počítačové simulace nezávislý na konkrétních hodnotách interakčních energií.Ty musí pouze odpovídat intuitivně chápaným relacím mezi nimi a musí nastavit systém pod kritickou teplotu fázového přechodu. Energie β i, j
Plyn x j = 0
Kapalina x j = 1
Vlákno x j = 2
Plyn xi = 0
− 40
− 10
Kapalina xi = 1
− 10
− 26
− 10
20
− 10
0
Vlákno xi = 2
20
Tab.1: V tabulce jsou uvedeny hodnoty interakčních konstant C ( xi , x j ) mezi uzly i a j , ve kterých mřížová proměnná X i a X j nabývá hodnot xi a x j . Při simulaci byly použity volné okrajové podmínky. Okolí každého uzlu sestávalo z krychle 3*3*3 uzly, ve kterém má uzel i 26 rovnocenných sousedů. Počáteční konfigurace byly tvořeny vláknem sestávajícím z řetízků uzlů seřazených těsně za sebou, který byl
Dílčí projekt: Systém projektování textilních struktur 3. Vývojové etapy rovnoběžný s osou x .Toto vlákno bylo obklopeno kapalinovým filmem neproměnného tvaru průřezu. Dynamika modelu byla typu Monte Carlo.Rovnoměrně náhodně byla vybrána dvojice uzlů a vypočítán hamiltonián původní konfigurace lišící se od původní právě výměnou zmíněných vybraných uzlů. Pro novou konfiguraci byl opět proveden výpočet hamiltoniánu. Nová konfigurace byla přijata tehdy, když její celková energie Q2 byla buď menší než celková energie konfigurace původní Q1 (Q2 〈Q1 ) , nebo když rovnoměrné náhodné číslo z ∈ 0,1 bylo menší než W1→2 = e
1 − ( Q2 −Q1 )
τ
. Veličina W1→2 je pravděpodobnost přechodu ze
stavu 1 do stavu 2. Pokud nebyly splněny výše popsané podmínky systém se vrátil do konfigurace původní. Dále pokračovala simulace opakováním právě popsaných kroků. Po provedení určitého počtu cyklů takových kroků byla výsledná konfigurace zaznamenána do paměti počítače. Za jednotku počtu provedených kroků byla zvolena MCSPS ( Monte Carlo Step per Site), což je takový počet provedených pokusů o výměnu dvojice uzlů, který je roven počtu uzlů v simulační krabici.V našem případě to je 54 *10 4 pokusů. Je známo [20], že výše popsaný algoritmus pro výměny částic přivede sledovaný systém do takové dráhy v konfiguračním prostoru, že pravděpodobnosti výskytu jeho konfigurací odpovídají svým rozložením pravděpodobnostem Pr rovnovážného systému G 1 Pr ( x ) = e Z
G −Q ( x )
τ
,
(26)
kde Z je statistická suma. Monte Carlo dynamika může nabývat různých podob co se týče vzdáleností, na které je povolena výměna mřížových proměnných. V našem případě jsme připouštěli výměnu dvou uzlů ( tak zvanou Kawasakiho dynamiku [20] ) jen na krátkou vzdálenost. Uzly navržené k výměně musely být sousedy. 3.6.7. VÝSTUPY POČÍTAČOVÉ SIMULACE Nejdůležitějšími výstupy počítačové simulace byly grafy, které znázorňovaly počty uzlů obsazené kapalinou v jednotlivých řezech simulační krabice, které byly kolmé na osu vlákna , tedy na osu x. Takto sestavené grafy citlivě odrážely systémy kapek, na které se vlivem Rayleighovy nestability rozpadaly původní stejnorodé kapalinové filmy. Ukázku řezu struktury modelu, řez je veden v ose vlákna , ukazuje obr.5. Časový vývoj počtu kapek je pak znázorněn na obr.6.
Obr.5: Na obrázku je znázorněno vlákno v řezu pokryté kapalinovým filmem, který je rozpadlý do jednotlivých kapek. Řez je veden v ose vlákna.
Dílčí projekt: Systém projektování textilních struktur 3. Vývojové etapy
Obr.6: Obrázek ukazuje časový vývoj tvorby kapek v MCSPS
Z obr.6 je patrné, že počet vytvořených kapek na vlákně z původního homogenního filmu kapaliny podléhá časovému vývoji. Počet kapek se postupně mění v důsledku slévání sousedních kapek.Vlnová délka λ není konstantní díky tomu, že na vlákně mezi kapkami existuje tenký film kapaliny prostřednictvím něhož si sousední kapky vyměňují kapalinové buňky. Tato situace je odlišná od rozpadu kapalinového sloupce čisté kapaliny, neobsahující vlákno, který vyústí v izolované ( nekomunikující) kapky. Časový vývoj počtu kapek středováný z pěti pokusů je zaznamenán na obr.7. Z výstupu simulace, které jsou graficky předvedeny na obr.7 vyplývá, že bylo možno jednoznačně určit počáteční počet kapek na vlákně v prvopočátcích úplného rozpadu kapalinového filmu. Výstupy počítačové simulace byly vzorkovány po 100 MCPS a byly zaznamenány konfigurace s pevně plně vyvinutými systémy kapek. Ranější konfigurace obsahovaly úseky nerozpadlého filmu a pozdější stádia vývoje jen plně vyvinuté a postupně slévající se kapky. Směrodatné odchylky v obr.7 nejsou tak vysoké v důsledku necitlivosti metody.
Dílčí projekt: Systém projektování textilních struktur 3. Vývojové etapy
Obr.7: Časový vývoj počtu kapek středovaný z pěti pokusů. Jednotkou času je MCSPS
Pro porovnání chování systému s vlnovou délkou λ Rayleighovy nestability (λ = 2.88πr ) jsme zvolili nejranější stadium vyvinutých kapek v počítačové simulaci. Mocnost kapalinového filmu na vlákně byla měněna. Jednotlivé průřezy počáteční konfigurace systémů s různými „poloměry“ kapalinového filmu ukazuje obr.8. Poloměr byl vypočítán tak, že ploše průřezu S sestávajícího se z určitého počtu uzlů ( S je počet uzlů násobený jednotkovou plochou v jednotkách (lattice unit) 2 = (l.u.) 2 byl přiřazen kruh o plošném obsahu πr 2 .
Obr.8 : Obrázek ukazuje příčné řezy vláken s počáteční konfigurací stejnorodého kapalinového filmu.Každému příčnému řezu je přiřazen poloměr kapalinového filmu S , kde S je plošný obsah řezu měřený v (l.u.) 2 poloměr r je vypočtený podle vztahu r =
π
vyjádřen v l.u.
Dílčí projekt: Systém projektování textilních struktur 3. Vývojové etapy Z nejranějších stádií rozpadu kapalinového filmu byly po uplynutí jistého počtu kroků simulace zjištěny počty vytvořených kapek N D . Rozdělením jejich počtu délkou vlákna L( L = 600l.u.) byla zjištěna vlnová délka Rayleighovy nestability příslušného systému, kterou dále budeme označovat λ sim . Taková vlnová délka byla porovnána s Rayleighovou vlnovou délkou λ , zjištěnou ze vzorce λ = 2.88πr . Každá simulace byla pro dané parametry opakována pětkrát. Výsledky těchto simulačních experimentů jsou přehledně zaznamenány v tab.2 a v obr.9 r [l.u ]
λ [l.u ] λ sim [l.u ]
2.034
2.585
2.821
3.038
3.432
3.785
18.405 19.667
23.413 23.393
24.180 25.523
27.489 26.818
31.050 31.383
34.243 32.329
Tab.2: Tabulka uvádí hodnoty vlnových délek λ a λ sim pro systémy s různým počátečním průměrem kapalinového filmu r .
Obr.8: Obrázek znázorňuje hodnoty Rayleighovy vlnové délky λ v závislosti na počátečním poloměru kapalinového filmu spolu s hodnotami vlnové délky λ sim zjištěnými pomocí použitého auto modelu.
Dílčí projekt: Systém projektování textilních struktur 3. Vývojové etapy 3.6.8. ZÁVĚR Ukázali jsme, že metoda Monte Carlo aplikovaná na vybraný Auto Model poskytuje výsledky, které jsou velmi blízké klasické teorii nestability kapalinových sloupců.V našem případě jsme sledovali kapalinové sloupce, které pokrývaly pevná rovná vlákna. Předpokládáme, že výsledky získané touto simulací by mohly být podobné s chováním reálných systémů, s nižším stupněm symetrie než je tomu u Rayleighovi nestability. Máme tím na mysli odchylky od osové symetrie systému způsobené vnějšími poli ( gravitačním nebo elektrickým) nebo odchylky způsobené nehomogenitou kapaliny při vyšetřování směsí nemísitelných tekutin jakými jsou například voda smíšena s oleji. Zdá se, že vyvinutá metoda by také mohla být vhodná pro studium chování kapalin na vláknech se složitou morfologií povrchu a chemickými nehomogenitami. Poděkování: Autoři tohoto příspěvku děkují Doc.RNDr. Aleši Linkovi CSc. za konzultace a inspiraci při formulování použitého modelu z pohledu teorie náhodných polí a Doc. Ing. Jiřímu Šedlbauerovi, Ph.D. za četné a cenné připomínky vedoucí k úpravě tohoto textu.
Literatura: 1. NEIMARK, A.V.: Thermodynamic equilibrium and stability of liquid films and droplets on fibers, J.Adhesion Sci. Technolog., Vol.13, No.10, pp.1137-1154 (1999). 2. RYONG-JOON ROE: Wetting of fine wires and fibres by a liquid film, Journal of Colloid and Interface Science, Vol. 50, No.1, pp.70-79. (1975). 3. PLATEAU, J.: Statique Experimentale et Theorique des Liquids Soumis aux Seules Forces Molecularics, Gauthier-Villars, paris (1873). 4. RAYLEIGH, Lord: Phil.mag. 34 , p.145 (1892). 5. TOMOTIKA, S.: Proc. Roy Soc. London Ser. A 153, p.302 (1936). 6. 7. MEISTER, B.J. and SCHEELE, G.F.: AICHE J. 13, p.682 (1967). 8. YOUNG, T.: In Miscellaneous Works, G.Peacock, ed. J.Murray, London, Vol. I, p.148 (1855) . 9. LAPLACE, P.S.: Mechanique Celeste, Supplement, Supplement to Book 10 (1806). 10. HILDEBRANDT,S. and TROMBA,A.: The shape and form in the natural world; The shape and form in the natural world; The parsimonious universe, Copernicus, SpringerVerlag, NewYork (1995). 11. BROCHARD, D.: Spreading of liquid drops on thin cylinder: The "manchon/droplet" transition, J.Chem.Phys., 84 , p.4664 (1986). 12. GRIGORJEV, A.I. -SHIR´AEVA, S.U.: Mechanism of elktrostatic polydispersion of liquid, J.Phys. D: Apply.Phys. 23 , pp. 1361-1370 (1990). 13. ADAMSON, A.W. and GAST, A.P.: Physical chemistry of surfaces, John Wiley & Sons, New York (1997).
Dílčí projekt: Systém projektování textilních struktur 3. Vývojové etapy 14. GEMAN, D.: Random fields and inverse problems in imaging, in Lecture, Notes in Mathematics, vol. 1427, pp. 113-193, Springer -Verlag (1991). 15. PAGET, R.D.: Nonparametric Markov random Filed Models for Natural texttile Images, Ph.D. thesis The University of Queensland, Australie, 16. BESAG, J.E.: Spatial interaction and the statistical of lattice systems, Journal of Royal Statistical Society, series B, vol. 36, pp. 192-326 (197). 17. MOUSSURIS,J.: Gibbs and Markov random systems with constraints, Journal of Statistical Physics, vol. 10, no.1, pp. 11-33 (1974). 18. HAMMERSLEY, J.M.- CLIFFORD,P.: Markov fields on finite graphs and lattices, unpublished (1971). 19. STAUFFER,D.: Scaling theory of percolation clusters, North-Holland publishing Co., Amsterdam (1979). 20. BINDER,K. and HEERMANN, D.W.: Monte Carlo Simulation in Statistical Physics, Springer-Verlag, Berlin (1997).
