Oblig 2: Dynamiske systemer og integrasjon

Læringsmål

I dissse oppgavene skal du lære og vise at du behersker følgende:

• Plotte og visualisere utviklingen i et dynamisk system.

• Simulere fartslover (ratelover).

• Sammenlikne eksperimentelle data med simuleringer.

• Forstå hva som menes med numerisk integrasjon, og gjøre rede for ulike metoder som kan brukes til dette.

Teoretisk bakgrunn

Vi kan også simulere kjemiske reaksjoner ved hjelp av modeller for reaksjonsfart. Disse modellene lar oss forutsi hvordan og hvor fort en kjemisk reaksjon vil gå. Dette kan brukes til å simulere alt fra industrielle prosesser til viktige reaksjoner i miljøet. Her skal vi se på modeller som kan forutsi hvordan det vil gå med ozonlaget i framtida.

Fartslover

Modeller for reaksjonsfart kaller vi fartslover. En fartslov beskriver endringen i konsentrasjon i en kjemisk reaksjon. La oss ta et enkelt eksempel der vi har to reaktanter og ett produkt:

\[A + B \rightarrow C\]

Fartsloven for denne reaksjonen må bestemmes eksperimentelt, og er derfor en empirisk lov. For eksempel kan endringen i konsentrasjonen til C være gitt ved:

\[\frac{d[C]}{dt} = k[A][B]\]

Her betyr \(\frac{d[C]}{dt}\) den deriverte av [C] med hensyn på tid (c’(t)). Det vil si at endringen i konsentrasjonen til produktet C er avhengig av konsentrasjonen til begge reaktanter i like stor grad. Men det kunne jo være at endringen i [C] varierte mer med [A] enn med [B], eller for eksempel ikke med [A] i det hele tatt. Da kunne vi henholdsvis fått disse modellene:

\[\frac{d[C]}{dt} = k[A]^2[B]\]

\[\frac{d[C]}{dt} = k[B]\]

Eksperimenter avgjør hvilken form vi gir fartslovene. Og dersom endringen av [C] er gitt ved \(\frac{d[C]}{dt} = k[A][B]\), kan vi ut fra reaksjonslikningen utlede følgende sammenhenger (forklar hvorfor!):

\[\frac{d[A]}{dt} = -k[A][B]\]

\[\frac{d[B]}{dt} = -k[A][B]\]

Fartslover for dannelse og nedbrytning av ozon i stratosfæren

Den såkalte Chapman-modellen kan benyttes for å simulere produksjon og nedbrytning av ozon i stratosfæren. Den er basert på følgende reaksjonslikninger med tilhørende reaksjonskoeffisienter:

\[O_2 \xrightarrow{uv} 2O\]

\[O_2 + O + M \rightarrow O_3 + M\]

\[O_3 \xrightarrow{uv'} O + O_2\]

\[O + O_3 \rightarrow 2O_2\]

hvor O, O\(_2\) og O\(_3\) er henholdsvis oksygen, dioksygen og ozon. M er en ikke-reagerende støtpartner, mens \(h \nu\) og \(h \nu '\) er energi tilført av UV-stråling med bølgelengde, \(\lambda\), under 242 nm og 336 nm, henholdsvis.

Den første reaksjonslikningen beskriver spaltingen av O\(_2\) til 2 O-atomer som resultat av UV-stråling. Den andre reaksjonslikningen viser den påfølgende reaksjonen mellom O\(_2\) og O som krever en kollisjon med M for å danne O\(_3\), mens de to siste reaksjonslikningene viser hvordan O\(_3\) brytes ned henholdsvis som resultat av UV-stråling for å danne O og O\(_2\), og gjennom reaksjon med O for produksjon av 2 O\(_2\)-molekyler.

Fartslovene for [O], [O\(_2\)] og [O\(_3\)] er gitt ved henholdsvis

\[\frac{d[\textrm{O}]}{dt} = 2 k_1 [\textrm{O}_2] - k_2 [\textrm{O}_2] [\textrm{O}] [\textrm{M}] + k_3 [\textrm{O}_3] - k_4 [\textrm{O}] [\textrm{O}_3\]

\[\frac{d[\textrm{O}_2]}{dt} = - k_1 [\textrm{O}_2] - k_2 [\textrm{O}_2] [\textrm{O}] [\textrm{M}] + k_3 [\textrm{O}_3] + 2 k_4 [\textrm{O}] [\textrm{O}_3]\]

\[\frac{d[\textrm{O}_3]}{dt} = k_2 [\textrm{O}_2] [\textrm{O}] [\textrm{M}] - k_3 [\textrm{O}_3] - k_4 [\textrm{O}] [\textrm{O}_3]\]

Ratekonstantene er gitt som følger:

\[k_1 = 3.0 \times 10^{-12} \textrm{ s}^{-1}\]

\[k_2 = 1.2 \times 10^{-33} \textrm{ cm}^6 \textrm{ molekyler}^{-2} \textrm{ s}^{-1}\]

\[k_3 = 5.5 \times 10^{-4} \textrm{ s}^{-1}\]

\[k_4 = 6.9 \times 10^{-16} \textrm{ cm}^3 \textrm{ molekyler}^{-1} \textrm{ s}^{-1}\]

Konsentrasjonene er gitt i molekyler per kubikkcentimerer (molec/cm\(^3\)). Steady-state-approksimasjonen sier at konsentrasjonen av intermediatene i en reaksjon er konstant, hvilket leder til uttrykkene nedenfor for [O\(_3\)] og [O] dersom vi regner med at [O\(_2\)] er konstant lik [O\(_2\)]\(_{t=0}\) under hele forløpet.

\[[\textrm{O}_3] = \sqrt{\frac{k_1 k_2}{k_3 k_4}}[\textrm{O}_2][\textrm{M}]^{\frac{1}{2}}\]

\[[\textrm{O}] = \frac{k_3[\textrm{O}_3]}{k_2[\textrm{O}_2][\textrm{M}]}\]

Ratekonstantene er gitt ved omtrent 25 km høyde, hvor \([\textrm{M}] \approx 9.0 \times 10^{17}\) molekyler cm\(^{-3}\). Systemet har følgende initialbetingelser:

\[[\textrm{O}_2]_{t=0} = 0.21[\textrm{M}]\]

\[[\textrm{O}]_{t=0} = 0\]

\[[\textrm{O}_3]_{t=0} = 0\]

2.1 Simulering av ozonnedbrytning

a) Lag et program som beregner og plotter [O\(_3\)] og [O] som funksjon av tid i intervallet \(t \in [0,100]\) ved å benytte Forward Euler-algoritmen på fartslovene i teoridelen med de gitte initialbetingelsene og tidssteg \(h = 0.001\). Plott med logaritmisk skala på \(y\)-aksen (plt.yscale(‘log’)).

b) Beregn og plott de samme verdiene med en backward-metode (‘BDF) ved å bruke funksjonen scipy.integrate.solve_ivp fra Scipy-biblioteket for \(t \in [0,10^8]\). Evaluer punktene for t = np.linspace(t0,tid_slutt,int(1E6)).

c) Sammenlikn de beregnede verdiene med BDF for [O\(_3\)] og [O] (de siste i arrayene) ved \(t=10^8\) s med verdiene fra steady-state-approksimasjonen. Beregn og kommenter det prosentvise avviket mellom steady-state og BDF-beregningene.

2.2 Numerisk integrasjon

Når vi løser fartslovene, integrerer vi dem. Det finnes ulike metoder for numerisk integrasjon, og hver metode har sine styrker og svakheter.

a) Forklar hvorfor fartslover er differensiallikninger. Hva finner vi når vi integrerer dem med Eulers metode?

b) Implementer rektangelmetoden og trapesmetoden for numerisk integrasjon som Python-funksjoner. Bruk begge metodene til å integrere en funksjon der du kjenner den analytisk integrerte. Sammenlikn de to metodene