04 · czerwiec 2026
Numeryczne modelowanie spalania metanu
Trzy własne solwery w Pythonie do spalania metanu w powietrzu, bez bibliotek kinetycznych: temperatura adiabatyczna (Newton–Raphson), zapłon w reaktorze 0D (RK4) i profil płomienia laminarnego 1D (różnice skończone, Newton i algorytm Thomasa).
Parametry
- Mieszanina
- CH₄ + powietrze, φ = 0,5 – 1,0
- Termodynamika
- wielomiany NASA-7
- Kinetyka
- 1-krokowa, Westbrook–Dryer
- Metody numeryczne
- Newton–Raphson · RK4 · MRS
- Temperatura adiabatyczna
- 2326 K (φ = 1)
- Pozorna energia aktywacji
- 175 kJ/mol
Cel projektu
Napisałem od zera trzy metody numeryczne do spalania metanu w powietrzu. Celowo nie użyłem gotowych bibliotek kinetycznych, takich jak Cantera. Bardziej niż sama fizyka spalania interesowała mnie strona numeryczna: dyskretyzacja, zbieżność i sprawdzenie, czy wynikom można ufać.
Globalna reakcja stechiometryczna:
CH₄ + 2 (O₂ + 3,76 N₂) → CO₂ + 2 H₂O + 7,52 N₂
Każda część to inny typ zadania numerycznego: równanie nieliniowe, układ równań różniczkowych zwyczajnych i równanie różniczkowe cząstkowe.
| Część | Zagadnienie | Metoda numeryczna |
|---|---|---|
| 1 | Temperatura adiabatyczna Tad(φ) | Newton–Raphson (równanie nieliniowe) |
| 2 | Zapłon w reaktorze 0D: T(t), Y(t) | Runge–Kutta 4. rzędu (układ ODE) |
| 3 | Profil płomienia laminarnego 1D | różnice skończone, Newton i algorytm Thomasa (PDE) |
Metodologia
Część 1: temperatura adiabatyczna
Temperaturę wyznaczam z bilansu entalpii mieszaniny metodą Newtona–Raphsona. Ciepło właściwe cp i entalpię h liczę z siedmiowspółczynnikowych wielomianów NASA w formacie CHEMKIN, osobno dla niskich i wysokich temperatur.
Dla mieszaniny stechiometrycznej metoda zbiega się w 3–4 iteracjach. Residuum bilansu entalpii spada z rzędu 10⁵ J poniżej 10⁻⁸ J, czyli kwadratowo, tak jak przewiduje teoria.
Część 2: reaktor 0D
Zapłon modeluję w jednorodnym, adiabatycznym reaktorze o stałym ciśnieniu, z jednokrokową globalną kinetyką Arrheniusa według Westbrooka i Dryera (1981). Rząd reakcji względem paliwa wynosi −0,3, a względem tlenu 1,3. Układ równań na YCH₄, YO₂ i T całkuję jawną metodą Rungego–Kutty 4. rzędu z krokiem 1 µs.
Czas zwłoki zapłonu policzyłem dla siedmiu temperatur początkowych z zakresu 1200–1500 K. Z wykresu Arrheniusa, ln(τ) od 1000/T, wyznaczyłem dopasowaniem liniowym pozorną energię aktywacji.
Część 3: płomień laminarny 1D
Równanie reakcji i dyfuzji rozwiązuję w układzie związanym z frontem płomienia:
ρ cp Su dT/dx = d/dx (λ dT/dx) + q ω(T)
Dyskretyzuję je metodą różnic skończonych na równomiernej siatce: człon dyfuzyjny schematem centralnym, a konwekcyjny przeciwprądowym (upwind). Powstały układ nieliniowy rozwiązuję iteracją Newtona–Raphsona. Jakobian jest trójdiagonalny, więc algorytm Thomasa robi każdy krok w O(N) operacji, a nie w O(N³), jak pełna eliminacja Gaussa.
Rozwiązanie sprawdziłem badaniem zbieżności siatkowej dla 51, 101, 201 i 401 węzłów.
Wyniki
| Wielkość | Wynik |
|---|---|
| Temperatura adiabatyczna, φ = 1 | 2326 K |
| Iteracje Newtona do zbieżności | 3–4 |
| Czas zwłoki zapłonu, T₀ = 1300 K | 18,3 ms |
| Pozorna energia aktywacji | 175 kJ/mol |
| Grubość frontu płomienia, N = 201 | 0,060 mm |
Pozorna energia aktywacji wyszła 175 kJ/mol, wyraźnie mniej niż 203 kJ/mol wpisane do modelu kinetycznego. Tego się spodziewałem: ujemny rząd reakcji względem paliwa spłaszcza zależność czasu zwłoki od temperatury.
Badanie siatkowe dało grubość frontu 0,090 / 0,070 / 0,060 / 0,055 mm dla 51, 101, 201 i 401 węzłów. Wartość maleje i zmienia się coraz wolniej, ale przy 401 węzłach jeszcze się nie ustaliła. Front płomienia mieści się na kilkunastu węzłach, więc wynik nadal zależy od gęstości siatki.
Ograniczenia modelu
Mechanizm jednokrokowy pomija dysocjację produktów, dlatego temperatura adiabatyczna wychodzi wyższa niż z pełnej równowagi chemicznej. Gdybym rozwijał projekt, zacząłbym od mechanizmu wielokrokowego. Potem dodałbym adaptacyjny krok całkowania dla układu sztywnego i wyznaczanie prędkości spalania laminarnego Su jako wartości własnej zagadnienia brzegowego.
Narzędzia
Python 3 (NumPy, Matplotlib) · LaTeX
Galeria wyników
