Mitteisotermilised reaktorid

Reaktsiooni kiirus on seotud energiatasemega reaktoris. Kui on rohkem energiat reaktoris, siis üldjuhul reaktsioon käib kiiremini. Reaktsioonid tekitavad või nõuavad samuti energiat ja võivad mõjutada, kui palju energiat oleks vaja anda või võtta reaktorist. Sellepärast, on mõnikord vaja mudelisse lisada energiabilansi.

Energia lisamine mudelisse

Energiabilansi on vaja arvesse võtta siis, kui temperatuur muutub või kui on vaja teada, kui palju energiat reaktor nõuab/toodab. Moolbilansist (ehk põhivõrrandist) me ei saa seda informatsiooni. Täpne energiabilanss sõltub ülesandest, aga üldiselt kuju on selline:

\[ \sum_i n_i C_{p,i} \frac{dT}{dt} = -\sum_i \dot n_{0,i} C_{p,i} (T - T_0) + r_A V \Delta H_{rxn} + \dot n_{faasimuutus} \Delta H_{faasimuutus} + Q \]

Kui koostada energiabilansi konkreetse ülesande jaoks, siis tuleb vaadata täpsemalt, mis energiaallikad ja energialäbikande liigid on olulised. Näiteks, kui on perioodiline reaktor ja on ainult üks faas, siis oluline on ainult reaktsiooni energia ja energialäbikanne reaktori ja ümbruse vahel.

Siin me tähistame energialäbikannet reaktori ja ümbruse vahel lihtsalt tähega Q, aga tegelikult võime kasutada soojuslevi valemeid selleks, et arvutada Q. Antud õppeaines me ei keskendu Q arvutamisele. Kuid see näitab, kuidas reaktsiooniprotsessid on seotud energia leviprotsessidega.

Energiabilanss torureaktoris

Kuna torureaktoris tingimused (k.a temperatuur) muutuvad mahuga, tasuks eraldi näidata torureaktori energiabilansi kuju. Torureaktori energiabilanssi on võimalik tuletada üldisest energiabilansist, aga selleks, et ei pea iga kord uuesti tuletama seda, näitame valemi üldkuju siin:

\[ \frac{dT}{dV} = \frac{ r_A \Delta H_{rxn} + \tilde n_{faasimuutus} \Delta H_{faasimuutus} + \tilde Q }{ \sum_i \dot n_i C_{p,i} } \]

Alustame üldise energiabilansiga. Kuna tavaliselt eeldame, et torureaktor töötab statsionaarses režiimis energiabilansi vasakupool on 0.

\[ 0 = \sum_i \dot n_{0,i} C_{p,i} (T - T_0) + r_A V \Delta H_{rxn} + \dot n_{faasimuutus} \Delta H_{faasimuutus} + Q \]

Nii nagu moolbilansiga, kui modelleerime torureaktorit tavaliselt tükeldame seda väikesteks mahtudeks ja käsitleme ühte tükki korraga. Ehk, ühes väikeses tükis meil on V asemel väike mahu muutus, dV. Samuti, erinevus siseneva voo (T₀) ja reaktori tüki temperatuuri vahel muutub väikeseks: dT. Pea meeles, et siin sisenev voog tuleb lihtsalt eelmisest tükist. Kaotame siis 0 moolkulu tähisest esimeses osas selleks, et oleks kergem meeles pidada, et see on moolkulu, mis siseneb teatud reaktori osasse. Nende muudatustega energiabilanss muutub selliseks:

\[ 0 = -\sum_i \dot n_i C_{p,i} dT + r_A dV \Delta H_{rxn} + \dot n_{faasimuutus} \Delta H_{faasimuutus} + Q \]

Kuna siin meid huvitab, kuidas temperatuur muutub mahuga reaktoris, jagame läbivalt dV-ga nii, et meil oleks \(\frac{dT}{dV}\), ehk temperatuuri tuletis mahu suhtes.

\[ 0 = -\sum_i \dot n_i C_{p,i} \frac{dT}{dV} + r_A \Delta H_{rxn} + \frac{\dot n_{faasimuutus}}{dV} \Delta H_{faasimuutus} + \frac{Q}{dV} \]

Pane tähele, et nüüd kõik osad näitavad energia voogu mahu kohta. Näiteks, oluline ei ole nüüd üldine soojusläbikande hulk (Q), vaid soojusläbikanne ühe mahu ühiku kohta. Asendame siis Q ja ̇n teiste tähistega aitamaks meil mäletada, et ühik on erinev.

\[ 0 = -\sum_i \dot n_i C_{p,i} \frac{dT}{dV} + r_A \Delta H_{rxn} + \tilde n_{faasimuutus} \Delta H_{faasimuutus} + \tilde Q \]

Nüüd saame teisendada valemit selleks, et liigutada temperatuuri tuletist vasakupoole.

\[ \sum_i \dot n_i C_{p,i} \frac{dT}{dV} = r_A \Delta H_{rxn} + \tilde n_{faasimuutus} \Delta H_{faasimuutus} + \tilde Q \]

\[ \frac{dT}{dV} = \frac{ r_A \Delta H_{rxn} + \tilde n_{faasimuutus} \Delta H_{faasimuutus} + \tilde Q }{ \sum_i \dot n_i C_{p,i} } \]

Siin valemi lugejas tähistame faasimuutuse moolkulu ja soojusläbikannet tildega kuna siin need väärtused on tegelikult mahu kohta. Kõik osad lugejas on mahu ühiku kohta kuna torureaktoris tahame teada, kuidas tingimused ja parameetrid muutuvad mahuga. Sellepärast, energiavood peavad olema mahu ühiku kohta.

Torureaktori energiabilansis, me lugejas sisuliselt liidame kokku kõik energiaallikad ja lahutame kõik energiakaod. Jagame seejärel süsteemi inertsiga, ehk kõikide ainete soojusmahtuvusega korda aine kogusega, selleks et arvutada temperatuuri muutumise kiirust.

Näide - patarei kuumenemine

Kui patareid on ebastabiilses seisundis, ühendid patareis võivad hakata reageerima, põhjustades tulekahju. Selline olukord võib näiteks tekkida, kui patarei temperatuur on liiga kõrge või kui laaditakse patarei liiga palju. On võimalik teha turvalisemaid patareisid, kui on mudelid, mis kirjeldavad neid protsesse patareis, mis võivad viia tulekahjuni.

Oletame, et liitium-ioon patarei on alguses toatemperatuuril ja siis pannakse keskkonda, mis on 240 kraadi juures. Kui kaua läheks enne, kui reaktsioonid patarei sees põhjustavad järsku temperatuuri tõusu (ja võibolla tulekahju)?

Sellises olukorras on mitu reaktsiooni, mis annavad soojust¹:

Metastabiilse tahke-elektrolüüdi liidese (ingl solid-electrolyte interface, SEI) lagunemine
Anoodi-elektrolüüdi reaktsioon
Elektrolüüdi lagunemine
Katoodi-elektrolüüdi reaktsioon

Eeldame, et patarei läbimõõt on 1,8 cm ja pikkus on 6,5 cm. Vajalikud kineetilised ja soojusläbikande parameetrid antakse Hewsoni¹ artiklis.

Alustame energiabilansiga. Meil on neli reaktsiooni, siis peame arvutama energiat, mida iga reaktsioon annab. Samuti toimub soojusläbikanne ümbrusega konvektsiooni ja kiirguse teel.

\[ \rho C_p \frac{dT}{dt} == Q_{rxn0} + Q_{rxn1} + Q_{rxn2} + Q_{rxn3} + Q_{konv} + Q_{kiirg} \]

Kõikide reaktsioonide jaoks kasutatakse sama vormi: reaktsiooni kiirust korda reaktsioonientalpiat korda aine mass mahu ühiku kohta (ρ). Kasutatakse aine massi kuna reaktsioonientalpiaid antakse massi ühiku kohta. Massid on antud mahu kohta kuna bilanss tervikuna antakse mahu ühiku kohta. Kuna kasutatakse erinevaid reaktsiooni kiirusevalemeid, näitame eraldi kõik neli valemit.

\[ Q_{rxn0} = k_0 Y_0 \rho_0 H_{rxn0} \] \[ Q_{rxn1} = k_1 Y_1 e^{-z/0.033} \rho_1 H_{rxn1} \] \[ Q_{rxn2} = k_2 Y_2 \rho_2 H_{rxn2} \] \[ Q_{rxn3} = k_3 Y_3 (1 - Y_3) \rho_3 H_{rxn3} \]

Soojusläbikanne konvektsiooni ja kiirguse teel võib arvutada järgmiste valemitega:

\[ Q_{konv} = h_{eff} \left( \frac{S}{V} \right) (T_{\infty} - T) \] \[ Q_{kiirg} = \epsilon \sigma \left( \frac{S}{V} \right) (T_{\infty} - T) \]

Siin h_eff on soojusläbikande tegur, \( \left( \frac{S}{V} \right) \) on pindala-mahu suhe ja T_∞ on ümbritseva keskkonna temperatuur.

Lisaks, peame koostama materjalbilanse iga aine jaoks. Kuna aineid ei lisata ega võeta patareist reaktsiooni ajal, patarei on sisuliselt perioodiline reaktor ja kontsentratsiooni muutus sõltub ainult reaktsiooni kiirusest. Siin kasutame normaliseeritud kontsentratsioon (Y), mis sisuliselt näitab, mis massosa antud ainest võib veel reageerida. Ehk, Y on umbes nagu konversiooniaste, aga on defineeritud massi järgi ja näitab, kui palju on järgi jäänud, mitte kui palju on reageerinud. Siin on massbilansid:

\[ \frac{dY_0}{dt} = -k_0 Y_0 \] \[ \frac{dY_1}{dt} = -k_1 Y_1 e^{-z/0.033} \] \[ \frac{dY_2}{dt} = -k_2 Y_2 \] \[ \frac{dY_3}{dt} = -k_3 Y_3 (1 - Y_3) \]

Nüüd võime kasutada Gekko paketti selleks, et lahendada diferentsiaalvõrrandit. Esiteks, paneme kirja reaktsiooni parameetrid, mida saime Hewsoni¹ artiklist. Metastabiilse tahke-elektrolüüdi liidese ρ väärtust ei olnud artiklis antud, siis lihtsalt eeldame, et on 100 kg m^-3.

Yalg_väärtused = np.asarray([0.15, 0.45, 1.0, 0.96]) # Y on sisuliselt normaliseeritud kontsentratsioon, ehk näitab, kui palju antud ainest on järgi jäänud
A_väärtused = np.asarray([1.67e15, 1.67e6, 5.1e25, 6.67e11]) # s^-1, sagedustegur
Ea_väärtused = np.asarray([135, 77.2, 274, 122]) * 1000 # J mol^-1, aktivatsioonienergia
Hrxn_väärtused = np.asarray([257.0, 1714, 155, 314]) * 1000 # J kg^-1, reaktsioonientalpia
rho_väärtused = np.asarray([100, 610.0, 1221, 407]) # kg m^-3, aine mass mahu ühiku kohta

Siis muudame neid Numpy massiive Gekko massiivideks. Saame seda teha for-tsükliga.

m = GEKKO()

ncomp = 4 # reaktsioonide arv
Yalg = m.Array(m.Param, (ncomp))
A = m.Array(m.Param, (ncomp))
Ea = m.Array(m.Param, (ncomp))
Hrxn = m.Array(m.Param, (ncomp))
rho = m.Array(m.Param, (ncomp))
for i in range(ncomp):
    Yalg[i].value = Yalg_väärtused[i]
    A[i].value = A_väärtused[i]
    Ea[i].value = Ea_väärtused[i]
    Hrxn[i].value = Hrxn_väärtused[i]
    rho[i].value = rho_väärtused[i]

Ja seejärel, defineerime ülejäänud parameetreid.

RGAS = m.Const(value=8.314462618) # J mol^-1 K^-1
rhoCp = m.Param(value=2.79e6) # J m^-3
pind_maht_suhe = m.Param(48.4) # m^-1
epsilon = m.Param(0.8)
sigma = m.Param(5.67e-8) # W m^-2 K^-4
heff = m.Param(value = 4.00377) # soojusläbikandetegur
Tvälis = m.Param(value=240) # K, ümbritseva keskkonna temperatuur

Siis peame Gekkole andma need ajahetked, mille juures sooviksime, et lahendust oleks arvutatud.

npts = 1000
t_lõpp = 3400
m.time = np.linspace(0, t_lõpp, npts) # s

Siis defineerime muutujaid, mida me tahame leida. Paneme ka vastavaid esialgseid väärtusi.

T = m.Var(25)
Y0 = m.Var(Yalg[0].value)
Y1 = m.Var(Yalg[1])
Y2 = m.Var(Yalg[2])
Y3 = m.Var(Yalg[3])

Nüüd saame kirja panna diferentsiaalvõrrandeid, mida me varem koostasime. Kuna terve energiabilanss on päris pikk, tükeldasime seda kasutades Gekko Intermediate objekte. Intermediate-iga on võimalik defineerida vahepealseid väärtusi või valemeid, mida hiljem kasutatakse lõplikus valemis (ehk Equation objektis).

z = m.Intermediate(0.033 + (Yalg[1] - Y1) + (Yalg[0] - Y0))
k = [None] * ncomp
for i in range(ncomp):
    k[i] = m.Intermediate(A[i] * m.exp(-Ea[i] / RGAS / (T + 273.15)))

Qrxn0 = m.Intermediate(k[0] * Y0 * rho[0] * Hrxn[0])
Qrxn1 = m.Intermediate(k[1] * Y1 * m.exp(-z/0.033) * rho[1] * Hrxn[1])
Qrxn2 = m.Intermediate(k[2] * Y2 * rho[2] * Hrxn[2])
Qrxn3 = m.Intermediate(k[3] * Y3 * (1 - Y3) * rho[3] * Hrxn[3])
Qkonv = m.Intermediate(heff * pind_maht_suhe * (Tvälis - T))
Qkiirg = m.Intermediate(epsilon * sigma * pind_maht_suhe * ((Tvälis + 273.15)**4 - (T + 273.15)**4))
m.Equation(T.dt() == (Qrxn0 + Qrxn1 + Qrxn2 + Qrxn3 + Qkonv + Qkiirg) / rhoCp)
m.Equation(Y0.dt() == -k[0] * Y0)
m.Equation(Y1.dt() == -k[1] * Y1 * m.exp(-z/0.033))
m.Equation(Y2.dt() == -k[2] * Y2)
m.Equation(Y3.dt() == -k[3] * Y3 * (1 - Y3))

Mudel on nüüd koostatud ja võime lahendada seda. Loome ka joonist selleks, et näha, kuidas temperatuur muutub ajas.

m.options.IMODE = 4
m.solve(disp=False)

plt.figure()
plt.plot(m.time, T.value)
plt.xlabel('Aeg (s)')
plt.ylabel("Temperatuur (°C)")
plt.savefig('aku-kuumenemine.png')
plt.show()

Joonisest võime näha, et 3200 sekundi juures (ehk umbes 53 minuti juures) patarei temperatuur hakkab järsult kasvama.

Teadlased ja insenerid kasutavad selliseid mudeleid selleks, et kiiresti hinnata, kas mõni uus aku keemia võiks olla piisavalt termiliselt stabiilne.

Pythoni fail, mida kasutasime siin lahendamiseks, on kättesaadav Githubist.

Viited

J. C. Hewson, “Understanding the limits of thermal runaway in lithium-ion battery systems,” Sandia National Lab. (SNL-NM), Albuquerque, NM (United States), SAND2016-6054C, Jun. 2016. Accessed: Oct. 21, 2021. [Online]. Available: https://www.osti.gov/biblio/1420927
H. S. Fogler, Essentials of Chemical Reaction Engineering. Pearson Education, 2010.