Vlastní implementace modelové funkce
Struktura programu v Matlabu bude stejná jako u ostatních programů: vždy potřebujeme hlavní program (skript) a modelovou funkci, která popisuje daný problém, tj. která z vektoru stavových proměnných vypočítá jejich derivace. Stavovými proměnnými budou v našem případě koncentrace složky A v jednotlivých KO (cA,1 až cA,N). Funkce bude obsahovat tyto části:
- rozbalení stavových proměnných a parametrů
- výpočet objemů KO a velikostí ploch mezi KO
- výpočet intenzit toků mezi jednotlivými KO
- bilance KO a z nich výpočet derivací stavových proměnných
Postupně teď ukážeme všechny části funkce pojmenované dfrpart.m a ty doplníme komentářem.
function [dydt] = dfrpart(~,y,pars) % Diffusion and reaction in spherical particle % ~ - first argument (time) is not needed by this function (autonomous ODE) % Y - vector of state variables (concentration of A) % DYDT - derivatives of state variables % PARS - model parameters % unpack model parameters csurf = pars(1); % surface concentration D = pars(2); % diffusion coefficient k = pars(3); % reaction rate constant R = pars(4); % particle radius N = max(size(y)); % get number of control volumes as the actual size % of the vector "y" |
Toto je klasický úvod pro soustavu ODE: Vstupními parametry je čas, vektor stavových proměnných a vektor parametrů. Parametry rozbalíme do proměnných s jasnými názvy. Poslední řádek je zajímavý. Počet KO odpovídá počtu stavových proměnných. Proto nepotřebujeme N zadávat zvlášť jako parametr, ale určíme jej z velikosti vektoru stavových proměnných y. Tento zápis bude fungovat správně ať už je vektor y na vstupu sloupcový nebo řádkový.
% calculate areas and volumes of control volumes r = zeros(N+1,1); Sr = zeros(N+1,1); V = zeros(N,1); h = R/N; % mesh size |
Prakticky jsme si vytvořili tři vektory: “r” bude vzdálenost hranic mezi KO od středu částice, “Sr” a “V” budou plochy mezi KO a objemy KO. Ze známeho poloměru částice a počtu KO určíme vzdálenost “h” mezi KO.
for i=1:N+1 r (i) = (i-1)*h; % radial position of the boundary Sr(i) = 4*pi*r(i)^2; % area between control volumes end for i=1:N V(i) = 4*pi/3 * (r(i+1)^3-r(i)^3); % volume of control volume end |
Ve dvou cyklech (for/end) jsme právě vypočítali prvky našich tří vektorů. Zápis i=1:N+1 znamená, že “i” nabývá hodnot od 1 do N+1. Všímněte si, že “r(1)” (a následně i “Sr(1)”) je rovno nule. Velikost plochy mezi KO vypočítáme jako plochu koule, objem KO vypočítáme jako rozdíl objemů dvou koulí. Protože jsme si vybrali, že “r(1)” bude nula, můžeme objem vnitřního objemu vypočítat stejným alogoritmem: od objemu koule odečítáme nulu.
% calculate the diffusive fluxes f = zeros(N+1,1); % mol/(s.m2) f(1) = 0.; % flux accros the center (not needed) f(2) = D / (1.5*h) * (y(1)-y(2)); for i=3:N f(i) = D / h * (y(i-1)-y(i)); end f(N+1) = D / (0.5*h) * (y(N)-csurf); % flux out of particle |
Tato část kódu počítá intenzitu toku. Protože uvažujeme jen difúzi, používáme Fickův zákon. Potřebujeme odhadnout derivaci koncentrace, pro kterou použijeme náhrady (Δc/Δr). Když se podíváme znovu pozorně na naše diskretizační schéma, vidíme, že pro toky f2 a fN+1 jsou jiné vzdálenosti (mezi zelenými křížky) než pro ostatní, tj. vnitřní toky. Tyto odlišnosti vyřešíme asymetrickými náhradami pro toky f2 a fN+1, u toků f3 až fN jsou náhrady symetrické.
Všimněte si, že toky budou kladné, když koncentrace blíže středu budou větší než koncetrace dále od středu: to je v souladu se směrem šipek, který jsme si určili.
% calculate concetration derivatives (control volume balances) dydt = zeros(N,1); % must be column vector for i=1:N dydt(i) = f(i)*Sr(i)/V(i) - f(i+1)*Sr(i+1)/V(i) - k*y(i); end end |
Poslední část kódu, kde byste měli identifikovat bilanci jednotlivých KO
Akumulace = Vstup – Výstup + Zdroj
a výpočet derivace koncetrace A podle času. Tyto derivace ukládáme do vektoru “dydt”, který musí být sloupcový. Dvě “end” na konci není chyba: první ukončuje cyklus, druhé celou funkci.
Ukázka hlavního programu
Teď je možné funkci vyzkoušet. V hlavním skriptu určíme parametry modelu, zadáme počáteční podmínky (tj. koncentraci v částici v čase nula) a zavoláme funkci pro integraci řešení soustavy ODE.
% diffusion and reaction in spherical particle (main script) cini = 1; % initial concentration inside particle csurf = 1; % concentration at the surface D = 1; % diffusion coefficient k = 10; % reaction rate constant R = 1; % particle radius N = 19; % number of control volumes pars = [csurf,D,k,R]; % pack model parameters c0 = ones(N,1) * cini; % initial condition tend = 0.2; % set end time of simulation [tt,yy] = ode15s(@(t,y) dfrpart(t,y,pars), 0:0.02:tend, c0); % integrate |
Výsledky simulace jsou uloženy v matici yy. Řádky matice odpovídají koncentračnímu profilu v jednotlivých časech, zatímco sloupce matice odpovídají časovému vývoji koncentrace v určitém místě částice (v určitém KO). Všimněte si, že jsme místo integračních mezí dodali přesné časy, ve kterých chceme znát hodnoty stavových proměnných ([0, 0.02, 0.04, ..., tend]) a to tak, aby mezi jednotlivými profily byl stejný časový krok. Protože se koncentrace mění jak s časem tak se vzdáleností od středu, je vizualizace výsledků komplikovanější než u obyčejných diferenciálních rovnic. Ukážeme proto jen jeden z příkladů, jak zakreslit vývoj koncentračního profilu v čase.
% plot concentration profile evolution subplot(1,2,1) h = R/N; rx = (1:N)*h - 0.5*h; rx(1)=0.0; plot([rx,R],[yy(1,:),csurf],'r-') hold on for i=2:max(size(tt))-1 plot([rx,R],[yy(i,:),csurf],'r-') end plot([rx,R],[yy(end,:),csurf],'bo-') hold off xlabel('radial coordinate') ylabel('concentration A') % plot concentration in the centre subplot(1,2,2) plot(tt,yy(:,1),'b--',tt,yy(:,end),'g-') xlabel('time') ylabel('concentration A') legend('centre','near surface') |
Na začátku je nutné do vektoru rx uložit pozice, ve kterých známe koncentrace (tj. vzdálenosti zelených křížků od středu částice). (1:N) je dvojtečkový operátor, který definuje vektor [1, 2, ..., N]. Ten se vynásobí vzdáleností mezi slupkami h a odečte se poloviční vzdálenost (tj. abychom se dostali do středu mezi hranicemi slupek). Pozici prvního elementu musíme separátně nastavit do středu částice. Při kreslení profilu v prvním časovém kroku zobrazíme i koncentraci na povrchu částice: na konec vektorů rx i yy(i,:) přidáme jeden prvek. Příkaz hold on nastaví chování příkazu plot tak, že před vykreslením čáry do grafu nesmaže již dříve nakreslené čáry. Poté v cyklu vykreslíme profily i ve zbývajících krocích. A nakonec popíšeme osy grafu.
A výsledek je zde:
Vidíme, jak se vyvíjí koncentrační profil (ubývá látky v částici) a jak vypadá koncentrační profil na konci simulace (modré kroužky). Na grafu vpravo pak pro kontrolu vidíme vývoj koncentrace v první (ve středu) a poslední slupce.
Print This Post
