Commit 87316f4ab8719050129f9f3ca31c21f68870b2b4

1D solver
Harjoitustyo/corrosive_darcy_flow.m
(408 / 0)
  
1% JM2011 HT (JY1).
2%
3% Tämä ohjelma ratkaisee yksinkertaisen 1+1D-mallin korrosiiviselle virtaukselle
4% huokoisessa väliaineessa. Väliaine oletetaan homogeeniseksi seokseksi, josta
5% osa syöpyy ja osa ei.
6%
7% vs = -(kappa / mu) p_x (Darcyn laki)
8% vs = phi v (efektiivinen virtausnopeus)
9%
10% rho_t + rho_x v + rho v_x = 0 (virtauksen kontinuiteetti)
11% rho = rho0 exp( beta [p - p0] ) (fluidin kompressoituvuus)
12%
13% cB_t = -k1 cB cA (konsentraatio B, huokoinen aine)
14% cA_t + vs cA_x = -k2 cB cA (konsentraatio A, syövyttävä, virtaava aine)
15%
16% phi = phi0 + ( (cB0 - cB(x,t)) / cB0 ) ( phi_infty - phi_0 ) (void fraction)
17%
18% missä
19%
20% vs = vs(x,t) = efektiivinen virtausnopeus (superficial fluid velocity)
21% v = v(x,t) = virtausnopeus tilanteessa ilman esteitä
22% kappa = huokoisen aineen permeabiliteetti, materiaalivakio
23% mu = virtaavan fluidin viskositeetti, vakio
24% p = p(x,t) = virtauksen painekenttä
25% rho = rho(x,t) = fluidin tiheys (ei vakio!)
26% beta = kompressoituvuus
27% rho0 = referenssitiheys kompressoituvuuden määrityksessä (huom. ei tarvita ratkaisussa!)
28% p0 = referenssipaine kompressoituvuuden määrityksessä (huom. ei tarvita ratkaisussa!)
29% cA = syövyttävän, virtaavan aineen konsentraatio
30% cB = huokoisen aineen konsentraatio
31% cB0 = huokoisen aineen konsentraatio alkutilanteessa
32% k1 = syöpymisreaktion reaktionopeus komponentille cB
33% k2 = syöpymisreaktion reaktionopeus komponentille cA
34% phi = phi(x,t) = tyhjän tilan osuus huokoisesta aineesta (void fraction)
35% 0 < phi < 1
36% phi_0 = tyhjän tilan osuus alussa
37% phi_infty = tyhjän tilan osuus, kun syöpyvä osa huokoisesta
38% aineesta on kulunut loppuun
39%
40% ja alaindeksit _t, _x merkitsevät osittaisderivaattoja.
41%
42% Näistä saadaan kolmen yhtälön kytketty systeemi (tarkemmin ks. raportti)
43%
44% p_t = g [p_x]^2 + 1/beta * [ g_x p_x + g p_xx ] (p_t)
45% cB_t = -k1 cB cA (cB_t)
46% cA_t = (kappa / mu) p_x * cA_x - k2 cB cA (cA_t)
47%
48% missä virtausvastus
49%
50% g(x,t) = (kappa / mu) * 1/phi(x,t)
51%
52% ja phi(x,t) lasketaan yllä annetulla kaavalla konsentraatiosta cB.
53%
54% Alkuehdot ovat
55%
56% cA(x,0) = 0
57% cB(x,0) = cB0
58% p(x,0) = p_infty - (phi_0 v_0 mu / kappa) * x
59%
60% ja reunaehdot
61%
62% p(0,t) = p_infty (paine vapaassa virtauksessa)
63% p_x(L,t) = 0 (tyypillinen ulosvirtausehto loppureunalla x = L)
64% cA(0,t) = cA0 (syövyttävän aineen vakiosyöttö alkureunalla x = 0)
65%
66% Näitä on eri määrä eri suureille, koska kertaluvut x:n suhteen ovat p: 2, cA: 1, cB: 0.
67%
68% missä p_infty on vakio (ympäristön paine). Paineen alkuehto seuraa
69% Darcyn laista integroimalla, kun oletetaan, että alkuhetkellä (toistaiseksi
70% ei-korrosiivinen) virtaus tapahtuu vakionopeudella v_0.
71%
72% Mallinnetussa tilanteessa siis systeemin läpi virtaa aluksi esim. pelkkää vettä
73% vakionopeudella (painegradientti = vakio), ja hetkellä t = 0 virtauksen
74% sekaan lisätty syövyttävä komponentti saavuttaa huokoisen aineen vasemman reunan.
75%
76% Ulosvirtausehdon takia joudutaan olettamaan, että tarkasteltava alue on "pitkähkö",
77% jotta virtaus ehtii tasaantua. Jotta voitaisiin tarkastella myös lyhyempiä "palikoita"
78% huokoista ainetta, täytyisi loppureunan reunaehto määritellä fysikaalisin perustein.
79%
80% JJ 2011-03-07
81function corrosive_darcy_flow
82 % Tehtäväparametrit tulevat tähän.
83 s = [];
84
85 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
86 % Parametrit (saa säätää vapaasti)
87 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
88
89 % Kvalitatiivinen demonstraatio. Hatusta vedetyt kertoimet.
90 % Voitaisiin käyttää myös fysikaalisesti järkeviä arvoja.
91 %
92 s.L = 10; % huokoisen "palikan" pituus, jota tarkastellaan
93
94 s.phi_0 = 0.5; % huokoinen aine = 50% tyhjää alussa
95 s.phi_infty = 0.9; % 1 - 0.9 = 10% syöpymätöntä ainetta
96 s.kappa = 1; % materiaalin permeabiliteetti
97
98 s.k1 = 0.2; % reaktionopeus, komponentti cB
99 s.k2 = 0.1; % reaktionopeus, komponentti cA
100 s.cB0 = 1; % huokoisen aineen syöpyvän osan konsentraatio alussa
101
102 s.mu = 1; % fluidin viskositeetti
103 s.beta = 1; % fluidin kompressoituvuuskerroin
104 s.p_infty = 20; % ympäristön paine
105 s.v_0 = 1; % vapaa virtausnopeus (ennen huokoiseen aineeseen törmäystä)
106
107 s.cA0 = 1; % syövyttävän aineen konsentraatio ennen törmäystä huokoiseen aineeseen
108
109 s.tmax = 10.0; % simulaation loppuaika
110 s.dt = 1e-3; % aika-askeleen koko samoissa yksiköissä kuin tmax
111 % (vakio vain yksinkertaisuuden vuoksi; RK4 ei rajoita)
112 s.nx = 51; % paikkadiskretoinnin solmujen lukumäärä
113
114 num_eps = 1e-8; % Numeerinen pyöristys, tätä pienemmät luvut = 0
115
116 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
117 % Laskentaskripti (särkyvää)
118 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
119
120 s.xx = transpose(linspace(0, s.L, s.nx)); % pystyvektori, helpompaa tämän koodin kannalta...
121 s.h = s.xx(2) - s.xx(1);
122
123 % Integroitavat suureet sijoitetaan yhteen pystyvektoriin.
124 % MATLAB-notaatiolla
125 %
126 % y = [ p ; cA ; cB ]
127 %
128 % missä jokainen suure on pystyvektori, jonka pituus on nx.
129
130 % Vapausasteita on siis yhteensä...
131 s.ny = 3*s.nx;
132
133 % ...ja suureiden indeksijoukot vektorissa y ovat
134 s.I_p = 1:s.nx;
135 s.I_cA = s.I_p + s.nx;
136 s.I_cB = s.I_p + 2*s.nx;
137 %
138 % Huomaa, että pisteet ovat järjestyksessä kasvavan x:n mukaan
139 % (tätä käytetään reunaehtojen käsittelyssä).
140
141 % Aika-askelten määrä (mukaanlukien alkuhetki t = 0)
142 s.nt = 1 + floor(s.tmax/s.dt);
143
144 % Varataan muisti tuloksille. NaN on kätevä arvo alustamattomalle osalle,
145 % koska MATLAB ei piirrä sitä kuviin.
146 %
147 % HUOM: aika pystyakselilla (rivi), suurevektorin y arvo vaaka-akselilla (sarake).
148 %
149 ys = nan(s.nt,s.ny);
150
151 % Ladotaan alkutila paikalleen:
152 %
153 y = nan(s.nx,1); % luodaan eksplisiittisesti pystyvektori (muutoin MATLAB arpoo jomman kumman)
154 y(s.I_p) = s.p_infty - (s.phi_0 * s.v_0 * s.mu / s.kappa) * s.xx;
155 y(s.I_cA) = 0;
156 y(s.I_cB) = s.cB0;
157 ys(1,:) = y'; % tulokset varastoidaan aika-askeleittain vaakavektoreina
158
159 fprintf(1, 'Ratkaistaan...\n');
160
161 % Integroidaan aika-askeleittain.
162 t = 0;
163 for k = 2:s.nt
164 % Reunaehdot paikalleen.
165 %
166 % Vasemman reunan Dirichlet-ehdot ovat helppoja.
167 % Pakotetaan arvot ennen integrointia:
168 y(s.I_p(1)) = s.p_infty;
169 y(s.I_cA(1)) = s.cA0;
170
171 % Oikealla reunalla joutuu hiukan temppuilemaan, koska tarvitaan p:n arvo,
172 % ja käytettävissä on vain p_x(L,t). Toisaalta, lisäksi käytettävissä
173 % on arvo edellisestä aika-askeleesta (tai 1. askeleella alkuarvo),
174 % ja yhtälöstä (p_t) saadaan p, kun tunnetaan p_x ja p_xx. Jälkimmäinen
175 % seuraa ensimmäisestä derivoimalla.
176 %
177 % Sijoitetaan siis näin saatava p_t(L,t) paikalleen funktioon f
178 % samalla tavalla kuin sisäosassakin, ja annetaan rk4:n hoitaa loput.
179 % Ks. funktion f toteutus.
180
181 % Aikaintegroidaan yksi aika-askel.
182 y = rk4(y, t, s.dt, @(y,t) f(s,y,t));
183
184 % Pyöristys, jottei numeerinen epätarkkuus kaada laskua.
185 % (Muutoin voi tulla hyvin pieniä negatiivisia arvoja konsentraatioille
186 % siellä, missä ne ovat lähellä nollaa.)
187 %
188 y( abs(y) < num_eps ) = 0;
189
190 % Tehdään pari sanity checkiä.
191 % Kaikki suureet p, cA, cB ovat 1) numeroita, 2) epänegatiivisia.
192 if sum(isnan(y)) || sum(y < 0)
193 % DEBUG
194 p = y(s.I_p)
195 cA = y(s.I_cA)
196 cB = y(s.I_cA)
197 k
198 t
199 error('kaboom');
200 end
201
202 % Tallennetaan tulos piirtoa varten.
203 %
204 % TODO: Tähän voisi virittää mekanismin, joka tallentaa vain joka n:nnen askeleen
205 % (tai yhteensä vain m kpl askelia, laskien automaattisesti n:n sopivasti).
206 %
207 ys(k,:) = y';
208
209 t = t + s.dt;
210 end
211
212 fprintf(1, 'Visualisoidaan...\n');
213
214 figure(1);
215 clf;
216
217 [X,T] = meshgrid(s.xx, 0:s.dt:s.tmax);
218
219 subplot(2,2, 1);
220% surf(X,T, ys(:,s.I_p), 'EdgeColor', 'none');
221% view(0,90);
222 contourf(X,T, ys(:,s.I_p));
223 xlabel('x');
224 ylabel('t');
225 title('p(x,t)');
226 colorbar('Location', 'SouthOutside');
227 colormap gray;
228 axis tight;
229
230 subplot(2,2, 2);
231 phi_ys = s.phi_0 + ( (s.cB0 - ys(:,s.I_cB)) / s.cB0 ) * ( s.phi_infty - s.phi_0 );
232 contourf(X,T, phi_ys);
233 xlabel('x');
234 ylabel('t');
235 title('\phi(x,t)');
236 colorbar('Location', 'SouthOutside');
237 axis tight;
238
239 subplot(2,2, 3);
240% surf(X,T, ys(:,s.I_cA), 'EdgeColor', 'none');
241% view(0,90);
242 contourf(X,T, ys(:,s.I_cA));
243 xlabel('x');
244 ylabel('t');
245 title('c_A(x,t)');
246 colorbar('Location', 'SouthOutside');
247 axis tight;
248
249 subplot(2,2, 4);
250% surf(X,T, ys(:,s.I_cB), 'EdgeColor', 'none');
251% view(0,90);
252 contourf(X,T, ys(:,s.I_cB));
253 xlabel('x');
254 ylabel('t');
255 title('c_B(x,t)');
256 colorbar('Location', 'SouthOutside');
257 axis tight;
258
259 fprintf(1, 'Valmis.\n');
260end
261
262% Yhtälön y_t = f(y, t) oikea puoli suurevektorille y.
263%
264% s = tehtäväparametrit (struct)
265% y = kentän arvo hetkellä t
266% t = aikamuuttujan arvo hetkellä t
267%
268function y_t = f(s, y, t)
269 % Void fraction phi(x,t). Huomaa pultattu laskentaruudukko suurevektorissa y.
270 % Eksplisiittistä riippuvuutta pelkästä aikamuuttujasta ei ole,
271 % joten aika t voidaan myös jättää pois.
272 %
273 phi = s.phi_0 + ( (s.cB0 - y(s.I_cB)) / s.cB0 ) * ( s.phi_infty - s.phi_0 );
274
275 % Virtausvastus g(x,t). Huomaa pultattu laskentaruudukko suurevektorissa y,
276 % joka implisiittisesti sisältyy suureeseen phi. Huomaa myös pultattu aika t.
277 %
278 g = (s.kappa / s.mu) * ( 1 ./ phi );
279
280 % Aikaevoluutioyhtälöt ovat siis
281 %
282 % p_t = g [p_x]^2 + 1/beta * [ g_x p_x + g p_xx ] (p_t)
283 % cA_t = (kappa / mu) p_x * cA_x - k2 cB cA (cA_t)
284 % cB_t = -k1 cB cA (cB_t)
285
286 % Lasketaan ensin apusuureita (kaikki, mitä tarvitaan enemmän kuin kerran).
287 p_x = ddx(y(s.I_p),s.h);
288
289 % Sisäosa:
290 y_t = [ g.*( p_x.^2 ) + 1/s.beta * ( ddx(g,s.h).*p_x + g.*d2dx2(y(s.I_p),s.h) ) ; ... % p_t
291 (s.kappa/s.mu) * p_x .* ddx(y(s.I_cA),s.h) - s.k2 * y(s.I_cB) .* y(s.I_cA) ; ... % cA_t
292 -s.k1*y(s.I_cA).*y(s.I_cB) ... % cB_t
293 ];
294
295 % Käsitellään paineen reunaehto alueen loppureunalla x = L.
296 %
297 % Tiedetään p_x(L,t) = 0. Siispä yhtälö (p_t) tulee muotoon
298 %
299 % p_t(L,t) = g/beta * p_xx(L,t)
300 %
301 % missä p_xx(L,t) = ddx(p_x(L,t)).
302 %
303 p_x(end) = 0; % reunaehdosta!
304 % (Tarvittavat pari sisempää pistettä seuraavat derivaattakentästä alueen sisällä, joka laskettiin jo yllä.)
305 p_xx_at_L_t = ddx_bwd(p_x((end-2):end),s.h); % tarvitaan vain loppupisteen arvo; ei lasketa turhaan alueen sisällä
306 p_xx_at_L_t = p_xx_at_L_t(end); % poimitaan laskun tulos vektorista
307
308 % Siispä reunaehto toteutuu, kun määräämme, että
309 y_t(s.I_p(end)) = (g(end)/s.beta) * p_xx_at_L_t;
310end
311
312% Standardi Runge-Kutta 4. Integroi yhtälöä
313%
314% y_t = f(y, t)
315%
316% yhden aika-askeleen. Suure y saa olla vektori, pysty- tai vaakaformaatilla
317% ei ole merkitystä.
318%
319% Parametrit:
320%
321% y0 = suureen y arvo aika-askeleen alussa
322% t0 = aika-askeleen alkuhetki
323% dt = aika-askeleen pituus
324% f = f(y, t) = yhtälön oikea puoli, function handle.
325%
326function y = rk4(y0, t0, dt, f)
327 k1 = f(y0, t0);
328 k2 = f(y0 + 1/2 * dt * k1, t0 + 1/2 * dt);
329 k3 = f(y0 + 1/2 * dt * k2, t0 + 1/2 * dt);
330 k4 = f(y0 + dt * k3, t0 + dt);
331 y = y0 + 1/6 * dt * (k1 + 2*k2 + 2*k3 + k4);
332end
333
334% Laskee skalaarikentän ensimmäisen paikkaderivaatan
335% toisen kertaluvun keskeisdifferenssillä tasavälisessä ruudukossa.
336% Päät käsitellään toisen kertaluvun forward- ja backward-differensseillä.
337%
338% Parametrit:
339% y = derivoitava kenttä (pystyvektoriformaatissa)
340% h = vektorin y pisteiden välimatka x-akselilla (vakio!)
341%
342function y_x = ddx(y, h)
343 % Keskiosa menee keskeisdifferenssillä vektoroidusti:
344 yp = [ y(2:end) ; NaN ];
345 ym = [ NaN ; y(1:end-1) ];
346 y_x = (yp - ym) / (2*h);
347
348 % Käsitellään päät toisen kertaluvun forward- ja backward-differensseillä.
349 % (Backward: http://www.geometrictools.com/Documentation/FiniteDifferences.pdf
350 % Forward samasta kaavasta etumerkit vaihtamalla.)
351 y_x(1) = (-3*y(1) + 4*y(2) - y(3)) / (2*h);
352 y_x(end) = ( 3*y(end) - 4*y(end-1) + y(end-2)) / (2*h);
353end
354
355% Laskee skalaarikentän ensimmäisen paikkaderivaatan
356% toisen kertaluvun takenevalla differenssillä tasavälisessä ruudukossa.
357%
358% Parametrit:
359% y = derivoitava kenttä (pystyvektoriformaatissa)
360% h = vektorin y pisteiden välimatka x-akselilla (vakio!)
361%
362% HUOM! Tämä funktio ei käsittele kahta ensimmäistä alkiota. Tässä ratkaisijassa
363% tätä käytetään vain paineen derivaattareunaehdon kanssa loppupäässä
364% ( p_xx(L,t) = d_dx( p_x(L,t) ), missä oikean puolen p_x saadaan suoraan reunaehdosta
365% ja p_xx lasketaan tällä funktiolla yhdistäen reunaehto sekä pari sisempää pistettä
366% alueen sisältä (joiden arvot tunnetaan edelliseltä aika-askeleelta) ).
367%
368function y = ddx_bwd(y, h)
369 ym1 = [ NaN ; y(1:end-1) ];
370 ym2 = [ NaN ; NaN; y(1:end-2) ];
371 y_x = ( 3*y - 4*ym1 + ym2 ) / (2*h);
372end
373
374% Laskee skalaarikentän toisen paikkaderivaatan
375% toisen kertaluvun keskeisdifferenssillä tasavälisessä ruudukossa.
376% Päät käsitellään toisen kertaluvun forward- ja backward-differensseillä.
377%
378% Parametrit:
379% y = derivoitava kenttä (pystyvektoriformaatissa)
380% h = vektorin y pisteiden välimatka x-akselilla (vakio!)
381%
382function y_xx = d2dx2(y, h)
383 % Sisäosa on jälleen helppo:
384 yp = [ y(2:end) ; NaN ];
385 ym = [ NaN ; y(1:end-1) ];
386 y_xx = ( ym - 2*y + yp ) / h^2;
387
388 % Toisen kertaluvun forward- ja backward-differenssit päille.
389 %
390 % Miten tämä tuotettiin:
391 % http://www.geometrictools.com/Documentation/FiniteDifferences.pdf
392 %
393 % Yhtälö (12) (vrt. example 1). d = 2, p = 2.
394 %
395 % Siispä toisen kertaluvun forward toiselle derivaatalle saadaan, kun
396 % A = [ 1 1 1 1 ; 0 1 2 3 ; 0 1 4 9 ; 0 1 8 27 ];
397 % b = [ 0 ; 0 ; 1 ; 0 ];
398 % C = A\b
399 % 2*C
400 % => C = 1/2 * (2, -5, 4, -1)
401 %
402 % Päiden käsittely tarvitaan tässä ratkaisijassa, koska muutoin
403 % NaNit propagoituisivat cB:hen (rk4:n sisällä, kun väliaskeleita lasketaan),
404 % ja saataisiin phi = NaN NaN ...
405 %
406 y_xx(1) = 0.5*y(1) - 2.5*y(2) + 2*y(3) - 0.5*y(4);
407 y_xx(end) = -0.5*y(end) + 2.5*y(end-1) - 2*y(end-2) + 0.5*y(end-3);
408end