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.

  • Bruke KI til å evaluere og forbedre egne programmer.

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\]

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) Bruk generativ KI til å forbedre programmet ditt, f.eks. ved å bruke funksjoner dersom du ikke har brukt det, eller ved å legge til relevant feilhåndtering. Kommenter kort hvilke endringer KI-en gjorde, og evaluer og begrunn om de er nyttige eller ikke.

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.

d) Forklar med ord hva som foregår i programmet ditt, både med utgangspunkt i programmeringen, matematikken og kjemien som danner grunnlag for simuleringen.