{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# II - Equilibrium calculations" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Equilibrium explanations" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Equilibrium calculations** are usually performed to obtain the adiabatic flame temperature, the equilibrium composition, and the thermodynamic state of a specific mixture under given conditions. These are virtually performed in every simulation.
\n", "For example, Cantera will call its equilibrium solver to initialize the gas state before trying to obtain a solution to the equations for a free flame. As such, it is interesting to understand how Cantera proceeds." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "There are 2 different types of solver currently implemented for equilibrium calculation in Cantera that\n", "deserves our attention :\n", "- **The element potential 'ChemEquil' solver**
\n", "The class ChemEquil implements a chemical equilibrium solver for single-phase solutions. It is a\n", "\"non-stoichiometric\" solver in the terminology of Smith and Missen, meaning that every intermediate\n", "state is a valid chemical equilibrium state, but does not necessarily satisfy the element constraints.\n", "
\n", "Non-stoichiometric methods are faster when they converge, but stoichiometric ones tend to be more\n", "robust.\n", "- **The 'VCS' chemical equilibrium solver**
\n", "The other type of solver is designed to be used to set a mixture containing one or more phases to\n", "a state of chemical equilibrium. It uses a \"stoichiometric\" algorithm, in which each intermediate\n", "state satisfies the element constraints but is not a state of chemical equilibrium. \n", "
\n", "More specifically, it\n", "implements the VCS algorithm, described in Smith and Missen, \"Chemical Reaction Equilibrium\". It\n", "finds a set of component species and a complete set of formation reactions for the non-components in\n", "terms of the components." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As expected, the **ChemEquil solver is the fastest** of the Cantera equilibrium solvers for many single-\n", "phase equilibrium problems (particularly if there are only a few elements but very many species), **but\n", "can be less stable**.
\n", "Problem situations include low temperatures where only a few species have non-zero\n", "mole fractions, precisely stoichiometric compositions (we will see an example shortly). In general, if\n", "speed is important, this solver should always be tried first before falling back to another one in case of\n", "failure." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### The function equilibrate" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The default setting in Cantera, when launching an equilibrium calculation without specifying\n", "the solver, is to try the 'vcs' before falling back to another vcs solver labelled 'gibbs' :\n", "\n", "gas.equilibrate('TP')\n", "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The equilibrate function can be applied on a single phase or on a mixture. Here, we recall its definition:\n", "\n", "equilibrate(self, XY, solver, double rtol, int maxsteps, int maxiter, int loglevel)\n", "\n", "
\n", "Parameters:\n", "\n", "XY - A two-letter string, which must be one of the set: ['TP','TV','HP','SP','SV','UV'].\n", "solver - Specifies the equilibrium solver to use. May be one of the following :\n", " element_potential = A fast solver using the element potential method.\n", " gibbs = A slower but more robust Gibbs minimization solver.\n", " vcs = The VCS non-ideal equilibrium solver.\n", " auto = The element potential solver will be tried first, then if it fails the gibbs solver will be\n", " tried.\n", "rtol - The relative error tolerance.\n", "maxsteps - Maximum number of steps in composition to take to find a converged solution.\n", "maxiter - This specifies the number of outer iterations on T or P when some property pair other than TP is\n", " specified (only for 'gibbs').\n", "loglevel - Determines the amount of output displayed during the solution process. 0 indicates no output, while larger numbers produce successively more verbose information.\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Useful python imports" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import cantera as ct\n", "import numpy as np\n", "import csv\n", "from matplotlib import *\n", "import matplotlib.pyplot as plt\n", "import sys" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1. Simple homogeneous equilibrium of a gas" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
Please create the gas object so that :
\n", "- the mechanism used is the gri30
\n", "- this is an air/CH4 mix
\n", "- the mix is at stoichiometry
\n", "- the temperature is 300K and the pressure 1 bar.\n", "" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "scrolled": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", " gri30:\n", "\n", " temperature 300 K\n", " pressure 1e+05 Pa\n", " density 1.1079 kg/m^3\n", " mean mol. weight 27.636 kg/kmol\n", " phase of matter gas\n", "\n", " 1 kg 1 kmol \n", " --------------- ---------------\n", " enthalpy -2.9979e+06 -8.2849e+07 J\n", " internal energy -3.0882e+06 -8.5344e+07 J\n", " entropy 7232.2 1.9987e+05 J/K\n", " Gibbs function -5.1676e+06 -1.4281e+08 J\n", " heat capacity c_p 1110.9 30700 J/K\n", " heat capacity c_v 810.01 22385 J/K\n", "\n", " mass frac. Y mole frac. X chem. pot. / RT\n", " --------------- --------------- ---------------\n", " H2O 0.12323 0.18904 -121.34\n", " CO2 0.15052 0.094518 -185.85\n", " N2 0.72625 0.71645 -23.38\n", " [ +50 minor] 2.887e-19 2.4934e-19 \n", "\n", "None\n" ] } ], "source": [ "gas = ct.Solution('gri30.yaml') # create an object representing the gas phase\n", "\n", " # set the composition of the gas\n", "\n", "gas.TP = 300,100000\n", "gas.set_equivalence_ratio(1.0,'CH4:1','O2:1, N2:3.79') # set initial state\n", "\n", "gas.equilibrate('TP') # equilibrate using Temperature (T) and Pressure (P)\n", "\n", "print(gas())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As you see, the gas has been equilibrated since it now shows only quantities for the product of the reaction (H2O and CO2). You can try to set yourself out of the perfect mixing (for example set CH4 to 0.4 and to 0.6) and see the impact on the species in the mix at the end." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2. Failure of a solver" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Cantera has 3 different equilibrium solvers, 2 of them are worth mentionning: \n", "- The 'ChemEquil' solver uses the element potential method for homogeneous equilibrium in gas mixtures. It is fast, but sometimes doesn't converge. \n", "- The 'VCS' solver uses the VCS algorithm (Gibbs minimization), which is slower but more robust. It can also handle multiple phases.
\n", "Here we'll solve a problem for which the ChemEquil solver fails, but the VCS solver has no problem." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Properties of the gas" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "pressure = 1.0e5 # pressure\n", "temperature = 400.0 # unburned gas temperature\n", "comp = 'CH4:0.5, O2:1, N2:3.76' # premixed gas composition\n", "\n", "gas = ct.Solution('gri30.yaml')\n", "gas.TPX = temperature, pressure, comp" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Initial state of the gas" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "******************************************************** \n", " Initial state :\n", "******************************************************** \n", "P = 1.0000e+05 Pa\n", "T = 4.0000e+02 K\n", "V = 1.2035e+00 m3/kg\n", "U = -2.6590e+05 J/kg\n", "H = -1.4554e+05 J/kg\n", "S = 7.5652e+03 J/kg/K\n", "\n", "\n" ] } ], "source": [ "print(\"******************************************************** \")\n", "print(\" Initial state :\")\n", "print(\"******************************************************** \")\n", "print(\"P = \", \"%10.4e \" % (gas.P) + \" Pa\")\n", "print(\"T = \", \"%10.4e \" % (gas.T) + \" K\")\n", "print(\"V = \", \"%10.4e \" % (gas.volume_mass) + \" m3/kg\")\n", "print(\"U = \", \"%10.4e \" % (gas.int_energy_mass) + \" J/kg\")\n", "print(\"H = \", \"%10.4e \" % (gas.enthalpy_mass) + \" J/kg\")\n", "print(\"S = \", \"%10.4e \" % (gas.entropy_mass) + \" J/kg/K\")\n", "print(\"\")\n", "print(\"\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Comparing chemical potentials and element potentials " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here, the chemical potentials (noted mu_xx) are compared to the corresponding calculated values with the element potentials (noted lambda_xx).\n", "For instance, mu_H2 = lambda_H x 2.\n", "This is a good way to check wether the solver has managed to compute the results correctly. The chemical potentials are the one of the vapor phase." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "Comparison between Chem potentials and element potentials:\n", "\n", "mu_H2 : -2.3501e+09 , -4.2513e+09\n", "mu_O2 : -8.8089e+07 , -4.2260e+09\n", "mu_OH : -2.3320e+09 , -4.2386e+09\n", "mu_H2O : -2.6153e+09 , -6.3643e+09\n" ] } ], "source": [ "chemeq = np.zeros(gas.n_species)\n", "chemeq = gas.chemical_potentials\n", "\n", "mu_H2 = chemeq[gas.species_index(\"H2\")]\n", "mu_OH = chemeq[gas.species_index(\"OH\")]\n", "mu_H2O = chemeq[gas.species_index(\"H2O\")]\n", "mu_O2 = chemeq[gas.species_index(\"O2\")]\n", "lambda_H = chemeq[gas.species_index(\"H\")]\n", "lambda_O = chemeq[gas.species_index(\"O\")]\n", "\n", "print()\n", "print(\"Comparison between Chem potentials and element potentials:\")\n", "print()\n", "s_mu_H2 = \"%11.4e\" % mu_H2\n", "s_lam_mu_H2 = \"%11.4e\" % (2.0 * lambda_H)\n", "print(\"mu_H2 : \", s_mu_H2, \", \", s_lam_mu_H2)\n", "\n", "s_mu_O2 = \"%11.4e\" % mu_O2\n", "s_lam_mu_O2 = \"%11.4e\" % (2.0 * lambda_O)\n", "print(\"mu_O2 : \", s_mu_O2, \", \", s_lam_mu_O2)\n", "\n", "s_mu_OH = \"%11.4e\" % mu_OH\n", "s_lam_mu_OH = \"%11.4e\" % (lambda_H + lambda_O)\n", "print(\"mu_OH : \", s_mu_OH, \", \", s_lam_mu_OH)\n", "\n", "s_mu_H2O = \"%11.4e\" % mu_H2O\n", "s_lam_mu_H2O = \"%11.4e\" % (2.0 * lambda_H + lambda_O)\n", "print(\"mu_H2O : \", s_mu_H2O, \", \", s_lam_mu_H2O)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Program equilibrate" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "By default, the method equilibrate can be used without giving the solver. It will try first the element_potential solver than the vcs solver. \n", "Just below, the default behavior of Cantera is illustrated. Don't use that in your everyday calculation, just call gas.equilibrate('TP'), Cantera will handle the rest." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "ChemEquil solver failed! Trying the vcs solver...\n" ] } ], "source": [ "try:\n", " # print(\"0\")\n", " gas.equilibrate('TP', solver='element_potential') # use the ChemEquil solver\n", "except:\n", " print(\"\")\n", " print(\"ChemEquil solver failed! Trying the vcs solver...\")\n", " gas.equilibrate('TP', solver='vcs', max_steps=1500)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Compare the results with the initial values" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "******************************************************** \n", " Final state :\n", " Tadiabatique = 400.0 K\n", "******************************************************** \n", "P = 1.0000e+05 Pa\n", "T = 4.0000e+02 K\n", "V = 1.2035e+00 m3/kg\n", "U = -3.0232e+06 J/kg\n", "H = -2.9029e+06 J/kg\n", "S = 7.5570e+03 J/kg/K\n", "\n", "\n" ] } ], "source": [ "print(\"\")\n", "print(\"******************************************************** \")\n", "print(\" Final state :\")\n", "print(\" Tadiabatique = \" + str(gas.T) + \" K\")\n", "print(\"******************************************************** \")\n", "print(\"P = \", \"%10.4e \" % (gas.P) + \" Pa\")\n", "print(\"T = \", \"%10.4e \" % (gas.T) + \" K\")\n", "print(\"V = \", \"%10.4e \" % (gas.volume_mass) + \" m3/kg\")\n", "print(\"U = \", \"%10.4e \" % (gas.int_energy_mass) + \" J/kg\")\n", "print(\"H = \", \"%10.4e \" % (gas.enthalpy_mass) + \" J/kg\")\n", "print(\"S = \", \"%10.4e \" % (gas.entropy_mass) + \" J/kg/K\")\n", "print(\"\")\n", "print(\"\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Comparing chemical and element equilibrium for the equilibrate mixture " ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "Comparison between Chem potentials and element potentials:\n", "\n", "mu_H2 : -1.7866e+08 , -1.7866e+08\n", "mu_O2 : -2.8958e+08 , -2.8958e+08\n", "mu_OH : -2.3412e+08 , -2.3412e+08\n", "mu_H2O : -3.2345e+08 , -3.2345e+08\n" ] } ], "source": [ "chemeq = gas.chemical_potentials\n", "mu_H2 = chemeq[gas.species_index(\"H2\")]\n", "mu_OH = chemeq[gas.species_index(\"OH\")]\n", "mu_H2O = chemeq[gas.species_index(\"H2O\")]\n", "mu_O2 = chemeq[gas.species_index(\"O2\")]\n", "lambda_H = chemeq[gas.species_index(\"H\")]\n", "lambda_O = chemeq[gas.species_index(\"O\")]\n", "\n", "print()\n", "print(\"Comparison between Chem potentials and element potentials:\")\n", "print()\n", "s_mu_H2 = \"%11.4e\" % mu_H2\n", "s_lam_mu_H2 = \"%11.4e\" % (2.0 * lambda_H)\n", "print(\"mu_H2 : \", s_mu_H2, \", \", s_lam_mu_H2)\n", "\n", "s_mu_O2 = \"%11.4e\" % mu_O2\n", "s_lam_mu_O2 = \"%11.4e\" % (2.0 * lambda_O)\n", "print(\"mu_O2 : \", s_mu_O2, \", \", s_lam_mu_O2)\n", "\n", "s_mu_OH = \"%11.4e\" % mu_OH\n", "s_lam_mu_OH = \"%11.4e\" % (lambda_H + lambda_O)\n", "print(\"mu_OH : \", s_mu_OH, \", \", s_lam_mu_OH)\n", "\n", "s_mu_H2O = \"%11.4e\" % mu_H2O\n", "s_lam_mu_H2O = \"%11.4e\" % (2.0 * lambda_H + lambda_O)\n", "print(\"mu_H2O : \", s_mu_H2O, \", \", s_lam_mu_H2O)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Saving the results" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Output written to 2-Output/all_mole_fractions.csv\n" ] } ], "source": [ "csv_file = '2-Output/all_mole_fractions.csv'\n", "with open(csv_file, 'w') as outfile:\n", " writer = csv.writer(outfile)\n", " writer.writerow(['phi', 'T (K)'] + gas.species_names)\n", " writer.writerow(['1', gas.T] + list(gas.X))\n", "print(('Output written to {0}'.format(csv_file)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 3. Perform adiabatic flame calculations as a function of equivalence ratio" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Set gas composition and interesting parameters" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will now perform **several constant pressure equilibrium calculations of a methane/air mixture, at\n", "an initial temperature of 300K and under atmospheric pressure**, in order to obtain the adiabatic flame temperature and burnt gas state for several equivalence ratios (from 0.3 to 10.0).
\n", "The goal of this exercise is to find a way to loop through several gaseous composition, in order to perform several computations in a single script and to plot the result." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Time to code a real cantera script !" ] }, { "cell_type": "markdown", "metadata": { "scrolled": true }, "source": [ "Then it is time for you to play ! Open another terminal folder and create your own Cantera python script.\n", "The goal is to compare the adiabatic flame temperature between gri 2.11 and gri3.0 (go search for them on https://chemistry.cerfacs.fr), for the following operating conditions :\n", "- a temperature of 400K\n", "- a pressure of 2 bars\n", "- a range of equivalence ratio between 0.3 to 3.5\n", "You should obtain the following results:\n" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGdCAYAAAAMm0nCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAA9hAAAPYQGoP6dpAABmzklEQVR4nO3deXxU9b3/8ddMMpmsk5VkiCQQ9kVkF+KCqCwCWrfeVqXqrba2FlqRe1FplYrWS1XqWqvXLtJeoa3+Wje0SNwAC7JHZd/CThJClsk+k5nz+2PISCBAEiZzsryfjwePOud8zzmf+TaaN2f5HIthGAYiIiIi7YjV7AJEREREmksBRkRERNodBRgRERFpdxRgREREpN1RgBEREZF2RwFGRERE2h0FGBEREWl3FGBERESk3Qk3u4DW4vP5OHLkCHFxcVgsFrPLERERkSYwDIPy8nLS09OxWs98nqXDBpgjR46QkZFhdhkiIiLSAgcPHqRbt25nXN9hA0xcXBzgnwCHw9HoGI/Hw7Jly5g4cSI2my2U5XVamvPQ05yHluY79DTnodeac+5yucjIyAj8Hj+TDhtg6i8bORyOswaY6OhoHA6HfuhDRHMeeprz0NJ8h57mPPRCMefnuv1DN/GKiIhIu6MAIyIiIu2OAoyIiIi0Ox32HhgREWkdhmFQV1eH1+s1uxTAfz9GeHg4NTU1baamju585jwsLIzw8PDzbnGiACMiIk3mdrs5evQoVVVVZpcSYBgGTqeTgwcPqu9XiJzvnEdHR9O1a1ciIiJaXIMCjIiINInP5yMvL4+wsDDS09OJiIhoE4HB5/NRUVFBbGzsWRufSfC0dM4Nw8DtdnPs2DHy8vLo06dPi/8/U4AREZEmcbvd+Hw+MjIyiI6ONrucAJ/Ph9vtJjIyUgEmRM5nzqOiorDZbOzfvz+wj5bQ/9MiItIsCglyvoLxM6SfQhEREWl3FGBERESk3VGAERERaecWLlxIQkKC2WWElAKMiIh0Cvn5+dx333307t2byMhI0tLSuPTSS3n55ZcDj4X36NEDi8WCxWIhOjqawYMH84c//KHBfj777DMsFgulpaVnPNbChQuxWCwMGDDgtHVvvvkmFouFHj16BPPrndPvf/97Lr/8chITE0lMTGT8+PGsXbv2rNscPXqU2267jb59+2K1Wpk5c+ZpY7Zs2cLNN98cmLvnnnuudb7AKRRg2ivDANcRXF++y75/Pkrh+nfMrkhEpM3au3cvw4YNY9myZfzP//wPmzZtYvXq1TzwwAMsWbKEjz76KDD2scce4+jRo2zevJnvfe97/PCHP+Rf//pXs48ZExNDYWEhq1evbrD8j3/8I5mZmef9nZrrs88+49Zbb+XTTz9l9erVZGRkMHHiRA4fPnzGbWpra+nSpQsPP/wwQ4YMaXRMVVUVPXv25Ne//jVOp7O1yj+NAkw7U7J7PXnPT8H1eA94ZgCOt26nx1fPkrjkLmrKi80uT0Q6GcMwqHLXhfyPYRjNqvMnP/kJ4eHhrF+/nu985zsMGDCAnj17cv311/P+++9z3XXXBcbGxcXhdDrp2bMnDz74IElJSeTk5DR7bsLDw7ntttv405/+FFh26NAhPvvsM2677bbTxr/zzjsMHz6cyMhIevbsybx586irqwusLy0t5Uc/+hFpaWlERkZy4YUXsmTJkgb7+PDDDxkwYACxsbFcc801HD16NLBu0aJF/OQnP2Ho0KH079+fP/zhD/h8Pj7++OMzfocePXrw/PPPc8cddxAfH9/omFGjRvH0009zyy23YLfbmzw/50t9YNqZw+/8kgvL/w2A17Cwy+hGqqWEJEsFu9a8x8Dxd5pcoYh0JtUeLwPnfhjy4259bBLREU37FXb8+PHAmZeYmJhGxzTWkM/n8/HWW29RUlLS4o6xd911F+PGjeP5558nOjqahQsXcs0115CWltZg3MqVK7njjjt44YUXuPzyy9mzZw/33HMPAL/85S/x+XxMnjyZ8vJyXn/9dXr16sXWrVsJCwsL7KOqqooFCxbwf//3f1itVr73ve/x3//93yxatKjR2qqqqvB4PCQlJbXou5lNZ2DambSKbQC8m/UIubdvodsvctmcci0AtduWmlmaiEibtHv3bgzDoF+/fg2Wp6SkEBsbS2xsLA8++GBg+YMPPkhsbCx2u51vf/vbJCYm8oMf/KBFxx42bBg9e/bk//2//4dhGCxcuJC77rrrtHHz5s3joYce4s4776Rnz55MmDCBxx9/nP/93/8F4KOPPmLt2rX885//ZMKECfTs2ZNrr72WyZMnB/bh8Xh45ZVXGDlyJMOHD2fGjBlnPbvy4IMPkp6ezvjx41v03cymMzDtSMXxw3QxjuMzLFx63fdJTkoGIHLgJFj5N7oXrwKfD9RkSkRCJMoWxtbHJply3PO1du1afD4f06ZNo7a2NrB89uzZ/Od//idHjx5l9uzZ/OQnP6F3795n3E9sbGzgn7/3ve/xyiuvNFh/11138dprr5GZmUllZSVTpkzht7/9bYMxX375Jf/+97954oknAsu8Xi81NTVUVVWRm5tLt27d6Nu37xnriI6OplevXoHPXbt2pbCwsNGxv/71r/nb3/7GZ5991uJOuGZTgGlHDm1ZTX9gv/UCsk6EF4ABoydSsSKSJEo5vH0NFwzMNq9IEelULBZLky/lmKV3795YLBZ27NjRYHnPnj0Bf2v7k6WkpNC7d2969+7Nm2++yeDBgxk5ciQDBw5sdP+5ubmBf3Y4HKetnzZtGg888ACPPvoot99+O+Hhp89XRUUF8+bN46abbjptXWRk5Gk1NsZmszX4bLFYGr1XaMGCBfz617/mo48+4qKLLjrnftsq/VW9HanIWwdAQUz/BsvjYmPZFjUcgPz174W8LhGRtiw5OZkJEybw29/+lsrKymZtm5GRwXe/+13mzJlzxjH1Yad3796kpqaetj4pKYlvfetbLF++vNHLRwDDhw9nx44dDfZV/8dqtXLRRRdx6NAhdu7c2az6T/XUU0/x+OOPs3TpUkaOHHle+zKbAkw7Yj/2NQCe1NMTc3X3qwFwHPo0pDWJiLQHv/vd76irq2PkyJH8/e9/Z9u2bezYsYPXX3+d7du3N7gZ9lT33Xcf7733HuvXr2/x8RcuXEhRURH9+/dvdP3cuXP5y1/+wrx589iyZQvbtm3jb3/7Gw8//DAAV1xxBWPHjuXmm28mJyeHvLw8/vWvf7F0adPvfXzyySd55JFH+NOf/kSPHj3Iz88nPz+fioqKwJg5c+Zwxx13NNguNzeX3NxcKioqOHbsGLm5uWzdujWw3u12B8a43W4OHz5Mbm4uu3fvbs4UNZsCTDvirNwOQGzW6ak5Y7T/EcBetduoKi0IaV0iIm1dr1692LRpE+PHj2fOnDkMGTKEkSNH8uKLL/Lf//3fPP7442fcduDAgUycOJG5c+e2+PhRUVEkJyefcf2kSZNYsmQJy5YtY9SoUYwZM4Znn32W7t27B8b84x//YNSoUdx6660MHDiQBx54AK/X2+QaXn75ZdxuN9/+9rfp2rVr4M+CBQsCY44ePcqBAwcabDds2DCGDRvGhg0bWLx4McOGDePaa68NrD9y5EhgzNGjR1mwYAHDhg1r8Y3PTWUxmvswfTvhcrmIj4+nrKys0WuS4L9j+4MPPmDKlCmnXTtsayqOHyb2xYH4DAsl9+0J3MBbzzAM9jw2hN7Gfr4e/RsGT27dH5yWak9z3lFozkOrI893TU0NeXl5ZGVltakbP30+Hy6XC4fDoTdlh8j5zvnZfpaa8vsbdAam3Ti0xd/Jcb/1gtPCC/hv1jra5XIA6nYuC2ltIiIioaYA005U7PNfez31Bt6TxVzo7wfQo2Q1hrfujONERETaOwWYdsJe+BXQ+A289QaMGo/LiCYRFwc2/ztUpYmIiIScAkw7cbYbeOtFRUWyI2YEAMc2LTnjOBERkfZOAaYdOLkDb+agMWcd6+7pbwmdcPizEFQmIiJiDgWYdqD+Bt4D1vRGb+A9WfeLrwegt2cn5cfP/Ip0ERGR9kwBph2ov4E3P2bAOcd2y8xip9X/Low9q95t1bpERETMogDTDjTlBt6TFaZeBoBv9yetVpOIiIiZmhVg5s+fz6hRo4iLiyM1NZUbbrihwcuxiouL+elPf0q/fv2IiooiMzOTn/3sZ5SVlTXYz4EDB5g6dSrR0dGkpqYye/Zs6uoaPvb72WefMXz4cOx2O71792bhwoUt/5btXFNu4D1ZTP8rAbjAtRE6Zp9CEZFWt2/fPiwWS4OXNUrb0awAs3z5cqZPn84XX3xBTk4OHo+HiRMnBl6OdeTIEY4cOcKCBQvYvHkzCxcuZOnSpdx9992BfXi9XqZOnYrb7WbVqlX8+c9/ZuHChQ1aNOfl5TF16lSuvPJKcnNzmTlzJj/4wQ/48MMPg/S124/m3MBbr9eIq3AbYaQZRRQe2HHuDURE5DQZGRkcPXqUCy+8EIDjx49zzTXXkJ6ejt1uJyMjgxkzZuByuRpsp7+Ah4hxHgoLCw3AWL58+RnHvPHGG0ZERITh8XgMwzCMDz74wLBarUZ+fn5gzMsvv2w4HA6jtrbWMAzDeOCBB4xBgwY12M93v/tdY9KkSU2urayszACMsrKyM45xu93G22+/bbjd7ibvN9S2LX/TMH7pMPIeHdCs7bY8PsYwfukwNrz1fCtV1jLtYc47Gs15aHXk+a6urja2bt1qVFdXm11KA16v1ygpKTG8Xm/Q9ln/++hkxcXFxu9+9ztj3bp1xr59+4yPPvrI6Nevn3HrrbcGxuzdu9eIjo42Zs2aZWzdutV48cUXjbCwMGPp0qVBq60tON85P9vPUlN+fxuGYYSfT/ipvzSUlJR01jEOh4PwcP+hVq9ezeDBg0lLSwuMmTRpEvfeey9btmxh2LBhrF69mvHjxzfYz6RJk5g5c+YZj1NbW0ttbW3gc30i9ng8eDyeRrepX36m9W1B+d51AORH9+eCZtRZ3GUUHN2KkbcSj+fe1iqv2drDnHc0mvPQ6sjz7fF4MAwDn8+Hz+czu5wA48Sl8vraGlNeXs69997LO++8g8PhYPbs2bz77rsMGTKEZ599lp49e3LXXXexa9cu3nnnHW688UZ++ctf0qtXLzZs2MDQoUOJj4/nRz/6UWCfGRkZ3HvvvSxYsCBw3JdffpmsrCyefvppAPr168fKlSt55plnmDBhQivPROg0Zc7PxufzYRgGHo/ntDeBN/XfnRYHGJ/Px8yZM7n00ksDp9dOVVRUxOOPP84999wTWJafn98gvACBz/n5+Wcd43K5qK6uJioq6rRjzZ8/n3nz5p22fNmyZURHR5/1u+Tk5Jx1vZm6HvQHmEM4KfjggyZvV+lzAnBB2QY+eP99sFhapb6Wastz3lFpzkOrI853eHg4TqeTiooK3G63f6FhQF21CcVEnfbftfLy8jMOv++++/j8889ZvHgxXbp0Yf78+WzcuJEBAwbgcrnw+XwsWLCABx54gOXLlwNQUVEBQGVl5WmXicD/5uY333yT7OzswPrPP/+cyy+/vMH4sWPHMmfOnEb30d6dbc7Pxu12U11dzYoVK067B7aqqqpJ+2hxgJk+fTqbN2/m888/b3S9y+Vi6tSpDBw4kEcffbSlh2myOXPmMGvWrAbHz8jIYOLEiWd9G3VOTg4TJkxos2+NLc29H4Csiydz0SXXNHm78vJLcT+/AKflOCMH9yY1s19rldgs7WHOOxrNeWh15Pmuqanh4MGDxMbGfvMGYXcl1l+fu8VDsPkeOgQRMYD/LEB5eTlxcXFYGvnLWnl5OX/96195/fXXue666wD4y1/+Qrdu3YiIiAi8Ufmqq67i5z//eWC7ffv2ARATE9Pg98htt93Gu+++S3V1Nddeey0LFy4MzEdRUREZGRkNxnfv3p3y8nJsNlujfwFvj8415+dSU1NDVFQUY8eObfRt1E3RogAzY8YMlixZwooVK+jWrdtp68vLy7nmmmuIi4vjrbfeavAvsdPpZO3atQ3GFxQUBNbV/2/9spPHOByOM/6fb7fbsdvtpy232Wzn/I9IU8aY4eQbeHsMvqRZNSYlpbDV1o+BdVvJ//pTLujV+Fkys7TVOe/INOeh1RHn2+v1YrFYsFqtWK0nngGxmtONw2q1Bo5dfwmjvrZT7du3D4/Hw5gxYwLrExMT6devX4NtRo0a1WD7+n9u8H2B5557jkcffZSdO3cyZ84c/vu//5vf/e53gfWn1nGm/bRn55rzc7FarVgslkb/PWnqvzfNCjCGYfDTn/6Ut956i88++4ysrKzTxrhcLiZNmoTdbufdd989LVllZ2fzxBNPUFhYSGpqKuA/1epwOBg4cGBgzAenXC7JyckhOzu7OeW2e0d2bqIvcNjqJOMcHXgbU5J6MRzZCvs+B34a9PpERLBFw8+PmHPcIIuJiWnSOKfTidPppH///iQlJXH55ZfzyCOP0LVr1xb9BVxaplmxafr06bz++ussXryYuLg48vPzyc/Pp7raf/3T5XIFHqv+4x//iMvlCozxer0ATJw4kYEDB3L77bfz5Zdf8uGHH/Lwww8zffr0wBmUH//4x+zdu5cHHniA7du387vf/Y433niD+++/P8hfv20rP7ITgCJ7Rou2j+k7DvDfB6N+MCLSKiwW/6WcUP9pxmWLnj17YrPZWLduXWBZWVkZO3fuPO+vX38mov4hkuzsbD7++OMGYzrjX8BDoVlnYF5++WUAxo0b12D5a6+9xn/+53+yceNG1qxZA0Dv3r0bjMnLy6NHjx6EhYWxZMkS7r33XrKzs4mJieHOO+/kscceC4zNysri/fff5/777+f555+nW7du/OEPf2DSpEkt+Y7tlu/4HgCqYru3aPtew6/E82kYTo5RcGAHad37B7M8EZF2IS4ujjvvvJPZs2eTlJREamoqv/zlLwOXMZrqgw8+oKCggFGjRhEbG8uWLVuYPXs2l156KT169AD8fwH/7W9/ywMPPMBdd93FJ598whtvvMH777/fSt+u82r2JaSzGTdu3DnHgP+GplMvETW2r02bNjWnvA4noizP/w9JPVu0fZwjgW22vgyo28ahTTkKMCLSaT3zzDP8+Mc/5tprr8XhcPDAAw9w8ODB025zOJuoqCh+//vfc//991NbW0tGRgY33XQTDz30UGCM/gIeOufVB0ZaV0LNIQCinH1bvI/iLhfD0W26D0ZEOrW4uDgWLVoU+FxZWcm8efMCbT7qnzg6WY8ePRr8pfzKK69k1apV5zyW/gIeGh3jdugOyPB5Savz3xiX1K3lj0DH9BsHnLgPRkSkk9q0aRN//etf2bNnDxs3bmTatGkAXH/99SZXJi2lANNGlebvJ9LiwWOE4eze8jMwvYZf5d+HcYyC/Xovkoh0XgsWLGDIkCGMHz+eyspKVq5cSUpKitllSQvpElIbdWz/VhKBo9ZUMhvpb9NUDe6DyV1GWve20dBORCSUhg0bxoYNOhPdkegMTBtVftT/eN/xFj5CfbLi1IsBMPIa75osIiLS3ijAtFG+Iv8j1NWxmee9r5i+VwC6D0ZERDoOBZg2KqJsHwBGUq/z3lf9fTBdjWMU7N9+3vsTkc6tKe0yRM4mGD9DCjBtVELNQQAi03qfY+S5xTkS2WPz3wh8aFPHe0OuiIRG/Ttqmvq2YJEzqf8ZOp/3hekm3jYo8Ai1BZIzgvOW1+OpF8ORbRh5K1E/GBFpibCwMBISEigsLAQgOjq6RW8iDjafz4fb7aampqbDvCyxrWvpnBuGQVVVFYWFhSQkJBAWFtbiGhRg2qDS/P0kBuER6pPF9h0HR/5MN9dG/3uR2sB/dESk/XE6nQCBENMWGIZBdXU1UVFRbSJQdQbnO+cJCQmBn6WWUoBpgwoPbAvKI9Qn6zXiatwn3ot0dP92uvYIzpkdEelcLBYLXbt2JTU1FY/HY3Y5AHg8HlasWMHYsWPP65KENN35zLnNZjuvMy/1FGDaoIoj/oZzxfZunP8zSH6xcfFsi+jHAM9WDm1cpgAjIuclLCwsKL+EgiEsLIy6ujoiIyMVYEKkLcy5Lha2QfWPULf0LdRnUpI6GgDLfvWDERGR9k0Bpg2qf4SaxJa9hfpMYvtdCUBG2QYMny+o+xYREQklBZg2KP7EI9R2Z5+g7rf3iCtxG2GkcZz8fduCum8REZFQUoBpYwyfF+eJt1AnZ/QP6r6jYxzsjvDf+3Jo07Kg7ltERCSUFGDamJKCAye9hTr4L14sO3EfjFX3wYiISDumANPGHNvvv7STb00lMkiPUJ8str//PphM10bdByMiIu2WAkwbU/8I9fGIbq2y/97Dr8RthNOFYg7v2dwqxxAREWltCjBtjLf+Eeq44D5CXS8qJpZdJ+6DOfql3oskIiLtkwJMGxPh2geAEeRHqE/mStN9MCIi0r4pwLQx8dX+R6ijgvAW6jOJG3gVAN3LdR+MiIi0TwowbcjJj1AnBfkR6pP1HjaOWsNGCqUc2v1Vqx1HRESktSjAtCElBQeJsripM6yt8gh1vcioGHbbdR+MiIi0Xwowbcix/VuBE49QR0a26rHK0sYAEKb7YEREpB1SgGlDyo/sBKColR6hPln8wKsByKrYiOHztvrxREREgkkBpg3xFe0GoDrIb6FuTO/hV1Bl2EnCxf5t61v9eCIiIsGkANOG1D9C7UtqvUeo69ntUeyOGgxA/pd6L5KIiLQvCjBtiKP6MABRqb1CcrzK9EsBiDz075AcT0REJFgUYNqQFG8BAPFdW68HTIPjXTQBgF6VudR53CE5poiISDAowLQRlWXHcVAJQEpGaAJMzwuzKSOGOEs1e7/SWRgREWk/FGDaiGOH/DfwFhNHfHxiSI4ZFh7Onuhh/uNuVj8YERFpPxRg2gjXUX+AKQpzhvS47szLAYg9siqkxxURETkfCjBtRE3RPgDKI7uG9LjOoRMB6F2zmdqaypAeW0REpKUUYNqKkgMA1Ma2fhO7k3XvO5RjJBJp8bB7w6chPbaIiEhLKcC0EREV/rdQWxIyQ3pci9XKfscIAMq3fRzSY4uIiLRUswLM/PnzGTVqFHFxcaSmpnLDDTewY8eOBmNqamqYPn06ycnJxMbGcvPNN1NQUNBgzIEDB5g6dSrR0dGkpqYye/Zs6urqGoz57LPPGD58OHa7nd69e7Nw4cKWfcN2wlFzFAB7l6yQH9vX4woAEgpWh/zYIiIiLdGsALN8+XKmT5/OF198QU5ODh6Ph4kTJ1JZ+c29E/fffz/vvfceb775JsuXL+fIkSPcdNNNgfVer5epU6fidrtZtWoVf/7zn1m4cCFz584NjMnLy2Pq1KlceeWV5ObmMnPmTH7wgx/w4YcfBuErt02BHjDO0DSxO1m34ZMA6O3eQYWrJOTHFxERaa7w5gxeunRpg88LFy4kNTWVDRs2MHbsWMrKyvjjH//I4sWLueqqqwB47bXXGDBgAF988QVjxoxh2bJlbN26lY8++oi0tDSGDh3K448/zoMPPsijjz5KREQEr7zyCllZWfzmN78BYMCAAXz++ec8++yzTJo0KUhfve0wowfMydJ79OOQxUk38tmzPochV30n5DWIiIg0R7MCzKnKysoASEpKAmDDhg14PB7Gjx8fGNO/f38yMzNZvXo1Y8aMYfXq1QwePJi0tLTAmEmTJnHvvfeyZcsWhg0bxurVqxvso37MzJkzz1hLbW0ttbW1gc8ulwsAj8eDx+NpdJv65WdaHyr5+7bTC38PmLjoWFPqOZwwkm4lS6ja/jGey29steO0lTnvTDTnoaX5Dj3Neei15pw3dZ8tDjA+n4+ZM2dy6aWXcuGFFwKQn59PREQECQkJDcampaWRn58fGHNyeKlfX7/ubGNcLhfV1dVERUWdVs/8+fOZN2/eacuXLVtGdHT0Wb9LTo65TdzchzfSCyigCys/+MCcGqyZjAaSC1fzQQhqMHvOOyPNeWhpvkNPcx56rTHnVVVVTRrX4gAzffp0Nm/ezOeff97SXQTVnDlzmDVrVuCzy+UiIyODiRMn4nA4Gt3G4/GQk5PDhAkTsNlsoSr1NBvf/AoKoTK6G1OmTDGlhuKCi+APv6Mv+0kcPYyE5NbpR9NW5rwz0ZyHluY79DTnodeac15/BeVcWhRgZsyYwZIlS1ixYgXdun3Tt8TpdOJ2uyktLW1wFqagoACn0xkYs3bt2gb7q39K6eQxpz65VFBQgMPhaPTsC4Ddbsdut5+23GaznXNymzKmNVnKDgHgjutmWh1p3bLYa+1BT98+Dm78iC5T7mrV45k9552R5jy0NN+hpzkPvdaY86bur1lPIRmGwYwZM3jrrbf45JNPyMpq+MjviBEjsNlsfPzxN/1EduzYwYEDB8jOzgYgOzubr7/+msLCwsCYnJwcHA4HAwcODIw5eR/1Y+r30dGY1QPmVIUpYwCo2/2JqXWIiIicS7MCzPTp03n99ddZvHgxcXFx5Ofnk5+fT3V1NQDx8fHcfffdzJo1i08//ZQNGzbw/e9/n+zsbMaM8f9ynDhxIgMHDuT222/nyy+/5MMPP+Thhx9m+vTpgTMoP/7xj9m7dy8PPPAA27dv53e/+x1vvPEG999/f5C/fttgZg+Yk0X199843a1kDRiGqbWIiIicTbMCzMsvv0xZWRnjxo2ja9eugT9///vfA2OeffZZrr32Wm6++WbGjh2L0+nkn//8Z2B9WFgYS5YsISwsjOzsbL73ve9xxx138NhjjwXGZGVl8f7775OTk8OQIUP4zW9+wx/+8IcO+Qg1mNsD5mS9R03EbYSRbhRyeO8WU2sRERE5m2bdA2M04W/lkZGRvPTSS7z00ktnHNO9e/dzPukybtw4Nm3a1Jzy2qWTe8B0yQx9D5iTxcTFs8U+kEHurzm04V9c0OtCU+sRERE5E70LyWTHDu0G/D1gHI5Ek6sBV9fLAIjYv9zkSkRERM5MAcZkrqP+AFMU5jS5Er/EwRMB6FW5Ee8p76cSERFpKxRgTFZTtA+A8sjW6bvSXL2HXI6LaBxUsuerttHjR0RE5FQKMGYrOQBAbWy3cwwMjXCbjT0xwwE4/tXSc4wWERExhwKMyQI9YBLN7QFzMnfmFQDEH9EZGBERaZsUYExW3wMmMsXcHjAnSx/hf51B79qtVFWUmluMiIhIIxRgTJbs9XckdnQ1twfMybr1HMhRuhBh8bJrnV6OJiIibY8CjIkqy4qJpwKALhnm9oA5mcVq5WDiaACqt39kcjUiIiKnU4Ax0bHDu4C20wPmZNbeVwGQdmy1yZWIiIicTgHGRK6je4C20wPmZL0unoLPsJDl209R/gGzyxEREWlAAcZENcfygLbTA+ZkiV26sje8JwD71r5vcjUiIiINKcCYyGhjPWBOdSz1EgCMPZ+aXImIiEhDCjAmslccAtpWD5iTxQ6cAECPsrUYPp/J1YiIiHxDAcZEcW2wB8zJ+owcT7URQRdKyNu6zuxyREREAhRgTJTiLQDaVg+Yk0VGxbAzaggAhZt0H4yIiLQdCjAmaas9YE5VnTkOgNhDK8wtRERE5CQKMCZpyz1gTtZ1xLUA9K35muoKl8nViIiI+CnAmKS+B8zxsDSTKzm7zD4XkU8XIix17Fynt1OLiEjboABjkpoi/yPULnvb6wFzMovVyoGkMQBUb9N7kUREpG1QgDGJr9T/CLUnpm0HGIDwvuMB6Fq0yuRKRERE/BRgTGKrPAKAz9E2m9idrNfoa6kzrHT3HSL/wC6zyxEREVGAMUtMtb8HjC05w+RKzi0+MYXdEf0BOLD2PZOrERERUYAxTULdMQDiunQ3uZKmKel6OQDheXqtgIiImE8BxgTeujpSfMcBSOja0+RqmibpomsA6F25njqP2+RqRESks1OAMUFxwQHCLT48RhgpaW3/HhiA3kPHUkosDqrYvWm52eWIiEgnpwBjguKjewEosiQRbrOZXE3ThIWHsyd2JAAlX6sfjIiImEsBxgSVhfsBKLG17SZ2p/L1vAqApKMrTa5EREQ6OwUYE3iKDwJQGek0uZLmybzY/1qB3p6dlB0vMLkaERHpzBRgzFDmb2JXF9v2m9idLK1bL/ZZMwmzGOxes8TsckREpBNTgDGBvcrfxM4S3/Z7wJwqv8ulAPh26rUCIiJiHgUYE8TW+i+/2NtBE7tTxQ6eAkDP0tX4vF6TqxERkc5KAcYESXWFAMSl9TC3kBboM2oClUYkyZSyd/Nqs8sREZFOSgEmxGqrK0jCBUByei+Tq2k+uz2KnTHDATi2SffBiIiIORRgQqzoSB4AVYadhKQuJlfTMu6sqwFIOKyGdiIiYg4FmBArPboPgGPWLlis7XP6u4++HoC+7m2UFR8zuRoREemM2udv0HasqsjfxM4VkWpyJS3nzOxDXv3j1KvfMbscERHphBRgQsxb4m9iVxXVvprYnaog1f92aj1OLSIiZmh2gFmxYgXXXXcd6enpWCwW3n777QbrKyoqmDFjBt26dSMqKoqBAwfyyiuvNBhTU1PD9OnTSU5OJjY2lptvvpmCgoadXQ8cOMDUqVOJjo4mNTWV2bNnU1dX1/xv2MZYyw8D4I1rHy9xPJPYC/1vp84q+0KPU4uISMg1O8BUVlYyZMgQXnrppUbXz5o1i6VLl/L666+zbds2Zs6cyYwZM3j33XcDY+6//37ee+893nzzTZYvX86RI0e46aabAuu9Xi9Tp07F7XazatUq/vznP7Nw4ULmzp3bgq/YtkRVHQUgLKF9B5i+oyZSaUSSosepRUTEBM0OMJMnT+ZXv/oVN954Y6PrV61axZ133sm4cePo0aMH99xzD0OGDGHt2rUAlJWV8cc//pFnnnmGq666ihEjRvDaa6+xatUqvvjiCwCWLVvG1q1bef311xk6dCiTJ0/m8ccf56WXXsLtdp/H1zWfw+0/0xSVkmlyJecnwh75zePUG98zuRoREelswoO9w0suuYR3332Xu+66i/T0dD777DN27tzJs88+C8CGDRvweDyMHz8+sE3//v3JzMxk9erVjBkzhtWrVzN48GDS0r55W/OkSZO499572bJlC8OGDTvtuLW1tdTW1gY+u1z+XisejwePx9NorfXLz7Q+6AyDFO8xsEBMSkbojttKanpcBVtXkXhkeZO/S8jnXDTnIab5Dj3Neei15pw3dZ9BDzAvvvgi99xzD926dSM8PByr1crvf/97xo4dC0B+fj4REREkJCQ02C4tLY38/PzAmJPDS/36+nWNmT9/PvPmzTtt+bJly4iOjj5rzTk5obkRta62kpstNQDkbtvL5t2HQ3Lc1lLrSSYb6OPezlv/eANbVGyTtw3VnMs3NOehpfkOPc156LXGnFdVVTVpXKsEmC+++IJ3332X7t27s2LFCqZPn056enqDsy7BNmfOHGbNmhX47HK5yMjIYOLEiTgcjka38Xg85OTkMGHCBGw2W6vVVm//tnWwFUqJ5VvXN34Jrr3Jm/8sWb4DdI+pZOg13znn+FDPuWjOQ03zHXqa89BrzTmvv4JyLkENMNXV1fz85z/nrbfeYurUqQBcdNFF5ObmsmDBAsaPH4/T6cTtdlNaWtrgLExBQQFOp//RYqfTGbhn5uT19esaY7fbsdvtpy232WznnNymjAmGqiL/I9RFYakkdJB/yQpSLycrfxHs/gib7Z4mbxeqOZdvaM5DS/Mdeprz0GuNOW/q/oLaB6b+fhPrKR1mw8LC8Pl8AIwYMQKbzcbHH38cWL9jxw4OHDhAdnY2ANnZ2Xz99dcUFhYGxuTk5OBwOBg4cGAwSw6pmuP+AFMRkXaOke1H7IWTAeipx6lFRCSEmn0GpqKigt27dwc+5+XlkZubS1JSEpmZmVxxxRXMnj2bqKgounfvzvLly/nLX/7CM888A0B8fDx33303s2bNIikpCYfDwU9/+lOys7MZM2YMABMnTmTgwIHcfvvtPPXUU+Tn5/Pwww8zffr0Rs+ytBe+Un+AqYlp303sTtZ31ATKc6JItpSx88uV9B0+zuySRESkE2h2gFm/fj1XXnll4HP9fSd33nknCxcu5G9/+xtz5sxh2rRpFBcX0717d5544gl+/OMfB7Z59tlnsVqt3HzzzdTW1jJp0iR+97vfBdaHhYWxZMkS7r33XrKzs4mJieHOO+/kscceO5/varrwCv9Nu7523sTuZBH2SDbHjWJ4xQqKN74LCjAiIhICzQ4w48aNwzCMM653Op289tprZ91HZGQkL7300hmb4QF0796dDz74oLnltWnR1f4nqCKSMkyuJLjqek2CL1eQcvQzs0sREZFOQu9CCqEEj/+enqiU7iZXElw9L7kBn2Ght3cPx47sM7scERHpBBRgQsTn9ZLiKwIgsWuWydUEV0paN3bZ+gKQt/otk6sREZHOQAEmRIqPHSbC4sVnWEjp2sPscoKu+AL/fVG2PctMrkRERDoDBZgQKTmyF4BjlkQi2vGTVGfSZfi3AOhXuYGa6kqTqxERkY5OASZEKgr3AVASnmpuIa2k1+BsCkki2lLLzjVLzS5HREQ6OAWYEHEfPwBAZWTHaWJ3MovVyr6kSwGo3vy+ydWIiEhHpwATIobL3wPGHZNuciWtxzbA35U34/jnGCc6L4uIiLQGBZgQiag8CoAlvuM0sTtV3zHXUmvYSDcKOLBjk9nliIhIB6YAEyKxNf4A09Ga2J0sJi6eHVFDATi6/h1zixERkQ5NASZEEur8PWBiUnuYW0grq+5xNQCOg5+YXImIiHRkCjAhUOdxk2yUAJDk7GFuMa0sY/SNAPSt3UJZ8TGTqxERkY5KASYEigsOEWYxqDOsJKVeYHY5rSo9qz951kzCLT52rVJXXhERaR0KMCFQWrAfgCJLEmHhzX5/ZruTnzbO/w871A9GRERahwJMCFQe8weY0vAUkysJjYRh/q68fcu/wOOuNbkaERHpiBRgQsBd4u8BU2XvmF14T9V3+JWU4MBBJTvW6d1IIiISfAowIWC4/I9Qu6OdJlcSGmHh4exO8HflrfjyPZOrERGRjkgBJgTCTzSxMxxdTa4kdMIGTAUg89hn6sorIiJBpwATAtE1BQDYEjr2E0gn63fJdYGuvPt3bDS7HBER6WAUYELA4fE3sYtKyTS5ktCJiUtge/QwAI6u1ePUIiISXAowrczw+Ujx+QNMfGrnCTAANVkTAUg8+LHJlYiISEejANPKykuKiLR4AEjp2t3kakKrxyU3AdDXs52i/IMmVyMiIh2JAkwrO56fB0AJcURGx5pcTWildevFrrDeWC0Geav+aXY5IiLSgSjAtLLywgMAlFiTTa7EHEUX+F/uGL5bXXlFRCR4FGBaWW3xIQDKI7qYXIk5UkfeAEC/yg3UVFWYW4yIiHQYCjCtrK7sCAA1UZ2jid2pel44hnxSiLbUsuOLJWaXIyIiHYQCTCsLq/AHGF9s52lidzKL1cr+lLEA1Gx+3+RqRESko1CAaWWRVf4mdmEJ6SZXYp6oC68FoGfxSnxer8nViIhIR6AA08pi3ccAsCdlmFyJefqNmUyFEUUXStjz5UqzyxERkQ5AAaaVJZ1oYhfXpfMGGHtkNDscYwAo3fSOydWIiEhHoADTimqqKkjA/+RNctce5hZjMqO//+WOFxR8anIlIiLSESjAtKLj+fsBqDYicCR2zseo6/W99CbcRhjdjUPUlB4xuxwREWnnFGBaUdmJAHPcmozF2rmn2pGQzPYo/8sdIws3mVyNiIi0d537t2orqzruf/9PmS3F5Erahuqe1wDQu3KDyZWIiEh7pwDTiupKDwNQHZlmciVtQ8/L/gOAQezm2JF95hYjIiLtmgJMK7KU++/18MR0zi68p+qS3oPt4f0B2L9aL3cUEZGWU4BpRbYTTewsjs7bxO5Uxd3GAxCT96HJlYiISHvW7ACzYsUKrrvuOtLT07FYLLz99tunjdm2bRvf+ta3iI+PJyYmhlGjRnHgwIHA+pqaGqZPn05ycjKxsbHcfPPNFBQUNNjHgQMHmDp1KtHR0aSmpjJ79mzq6uqa/w1NFFNTCEBEYjeTK2k7nBffBED/mlxcpcdNrkZERNqrZgeYyspKhgwZwksvvdTo+j179nDZZZfRv39/PvvsM7766iseeeQRIiMjA2Puv/9+3nvvPd58802WL1/OkSNHuOmmmwLrvV4vU6dOxe12s2rVKv785z+zcOFC5s6d24KvaJ6EOn8X3piUztvE7lQZfS4ij3QiLF52/luXkUREpGXCm7vB5MmTmTx58hnX/+IXv2DKlCk89dRTgWW9evUK/HNZWRl//OMfWbx4MVdddRUAr732GgMGDOCLL75gzJgxLFu2jK1bt/LRRx+RlpbG0KFDefzxx3nwwQd59NFHiYiIaG7ZIeetqyPZKAELJDi7m11Om7ItcgRZNUewbH8fpv7Q7HJERKQdanaAORufz8f777/PAw88wKRJk9i0aRNZWVnMmTOHG264AYANGzbg8XgYP358YLv+/fuTmZnJ6tWrGTNmDKtXr2bw4MGkpX3z9M6kSZO499572bJlC8OGDTvt2LW1tdTW1gY+u1wuADweDx6Pp9F665efaf35KMo/QFeLD69hwZHkbJVjtEcej4fylOFw6D36ub6gotyFPTLK7LI6tNb8OZfTab5DT3Meeq05503dZ1ADTGFhIRUVFfz617/mV7/6FU8++SRLly7lpptu4tNPP+WKK64gPz+fiIgIEhISGmyblpZGfn4+APn5+Q3CS/36+nWNmT9/PvPmzTtt+bJly4iOjj5r3Tk5OU39ik1WU7SX7wLHSWB1zkdB3397FpGcReHBRFItJfz9/54n8oKLzC6pU2iNn3M5M8136GnOQ6815ryqqqpJ44J+Bgbg+uuv5/777wdg6NChrFq1ildeeYUrrrgimIdrYM6cOcyaNSvw2eVykZGRwcSJE3E4HI1u4/F4yMnJYcKECdhstqDW89VHi+EglNi6MGXKlKDuuz2rn/M9yWNJLX6HHp7tDJ/ykNlldWit+XMup9N8h57mPPRac87rr6CcS1ADTEpKCuHh4QwcOLDB8gEDBvD5558D4HQ6cbvdlJaWNjgLU1BQgNPpDIxZu3Ztg33UP6VUP+ZUdrsdu91+2nKbzXbOyW3KmOaqcx0FoDIiVf9CNSL6ohvgs3foU7ICq8VCWHhQfxSlEa3xcy5npvkOPc156LXGnDd1f0HtAxMREcGoUaPYsWNHg+U7d+6ke3f/jawjRozAZrPx8ccfB9bv2LGDAwcOkJ2dDUB2djZff/01hYWFgTE5OTk4HI7TwlFbZZT5u/B6otWFtzF9Rk2kjBiScLF97TKzyxERkXam2X/traioYPfu3YHPeXl55ObmkpSURGZmJrNnz+a73/0uY8eO5corr2Tp0qW89957fPbZZwDEx8dz9913M2vWLJKSknA4HPz0pz8lOzubMWPGADBx4kQGDhzI7bffzlNPPUV+fj4PP/ww06dPb/QsS1sUVum/V8enJnaNskXY2ZIwllGl/6J80z/hEl1mExGRpmv2GZj169czbNiwwJNAs2bNYtiwYYEeLTfeeCOvvPIKTz31FIMHD+YPf/gD//jHP7jssssC+3j22We59tprufnmmxk7dixOp5N//vObniBhYWEsWbKEsLAwsrOz+d73vscdd9zBY489dr7fN2SiTjSxC4+/wORK2i7bhdcDkHXsE3xer8nViIhIe9LsMzDjxo3DMIyzjrnrrru46667zrg+MjKSl1566YzN8AC6d+/OBx980Nzy2ox4tz/ARKWoC++Z9L/0W1SunEma5Tg7v1xJ3+HjzC5JRETaCb0LqRUYPh/JPn+bfEeqmtidSWRUDNsd/vuejq970+RqRESkPVGAaQXlrhKiLf6meinpPcwtpq0b4L+MlJH/McaJx/BFRETORQGmFZQczQPARQzRMY33oBG//pffSI1ho5txlLyt68wuR0RE2gkFmFbgKvS/efu4NcXkStq+mLgEtsWMAqBgjS4jiYhI0yjAtILq4kMAlEd0MbmS9qGu77UAOA+rH4yIiDSNAkwr8JYdAaAmMtXkStqHvmO/g8cII8u3n4O7vjS7HBERaQcUYFqBtcL/2gNvjLrwNkV8Uhe2RQ0F4NCqN8wtRkRE2gUFmFYQUe3vAWN1dDW5kvajutdUAJIP6jKSiIicmwJMK4iuPQaALUGvEWiqXpd/B59hoW/dTo7u33HuDUREpFNTgGkF8XXFAEQn6zUCTZXizGCb/UIA9n/+N5OrERGRtk4BJsgMn48kwx9gElIzTK6mfano5X8aKTGv/b5CQkREQkMBJshKjxcQYfG/mDApTQGmOXpdfis+w0K/uu3kH9x97g1ERKTTUoAJspICfxO7EhxE2CNNrqZ9SUnvzvaIQQDsW7HY5GpERKQtU4AJsooifxO7EmuSyZW0T66e/qeREvb9y+RKRESkLVOACbKaksMAVEToNQIt0XPsrQD092yl4NAek6sREZG2SgEmyHxl+YC68LZU6gVZbLMNBCBvxV9NrkZERNoqBZggs1QcBcAbowDTUmVZ/stI8Xnvm1yJiIi0VQowQWar9jexs8SpC29LZZ24jDTAs5XCw3kmVyMiIm2RAkyQxZzowmtPVIBpqbRuvdgePgCAvbqMJCIijVCACTJH3XEAopK6mVxJ+1aaNQUAx15dRhIRkdMpwASR4fORbJQAEJ+qAHM+eoy9DYD+7i0UHdlvcjUiItLWKMAEUVlxIRGWOkBdeM+XM6M3O8L7Y7UY7FFTOxEROYUCTBCVFPjPFJQQhz0y2uRq2r+SHpMBiNuzxORKRESkrVGACaKKIn8TuxJrssmVdAw9rvge4L+MpKeRRETkZAowQVRdXN+FVwEmGJwZvdluG4jVYrB3+etmlyMiIm2IAkwQeV3+Jna16sIbNKU9rwUgca8uI4mIyDcUYILIWu5/jUBdtAJMsPS6Yho+w0K/uu0c3b/D7HJERKSNUIAJovouvFaHmtgFS5f0HmyzXwjA/hWLTK5GRETaCgWYIIo+0YXXlqAAE0wVva8HIHmfmtqJiIifAkwQObz+LrzRyWpiF0x9xt1GnWGlj3c3h3ZvNrscERFpAxRggsTw+UjxFQPg6KImdsGUlHoB2yKHAnDoczW1ExERBZigcZUcC3ThTXYqwARbdV//ZaTUA7qMJCIiCjBBU1xwAIBSYtWFtxX0G3crHiOMnr597N++0exyRETEZAowQVJRdAhQF97WEp+cxtboEQAc+bcuI4mIdHYKMEFSU9+F16YA01rc/W8EIP3wvzB8PpOrERERMynABEldmb8Lb4268Laa/uNuodaw0d13iLyt68wuR0RETNTsALNixQquu+460tPTsVgsvP3222cc++Mf/xiLxcJzzz3XYHlxcTHTpk3D4XCQkJDA3XffTUVFRYMxX331FZdffjmRkZFkZGTw1FNPNbfUkLJWnOjCG6MA01ri4pPYGjsagIJVejeSiEhn1uwAU1lZyZAhQ3jppZfOOu6tt97iiy++ID09/bR106ZNY8uWLeTk5LBkyRJWrFjBPffcE1jvcrmYOHEi3bt3Z8OGDTz99NM8+uijvPrqq80tN2Rs1YUAWOLUxK41+QZ9G4AeR3QZSUSkMwtv7gaTJ09m8uTJZx1z+PBhfvrTn/Lhhx8yderUBuu2bdvG0qVLWbduHSNHjgTgxRdfZMqUKSxYsID09HQWLVqE2+3mT3/6ExEREQwaNIjc3FyeeeaZBkGnLYmuLQIgQl14W9Wgcf9B5ZqH6Go5xvb1H9P/4glmlyQiIiYI+j0wPp+P22+/ndmzZzNo0KDT1q9evZqEhIRAeAEYP348VquVNWvWBMaMHTuWiIiIwJhJkyaxY8cOSkpKgl1yUDjq/AFGXXhbV2R0LFsTxgJQtlZPI4mIdFbNPgNzLk8++STh4eH87Gc/a3R9fn4+qakN7xMJDw8nKSmJ/Pz8wJisrKwGY9LS0gLrEhMTT9tvbW0ttbW1gc8ulwsAj8eDx+NptJb65Wda31SGz0eyrwQsEJ3kPO/9dWTBmPPwi/4DVi6jT9FHVFdVEm6LOPdGnViwfs6laTTfoac5D73WnPOm7jOoAWbDhg08//zzbNy4EYvFEsxdn9P8+fOZN2/eacuXLVtGdPTZG8vl5OSc17Hraiq42eKf8PVfbidsy97z2l9ncD5z7vOG0d2II8ni4u9/fpbI9MFBrKzjOt+fc2kezXfoac5DrzXmvKqqqknjghpgVq5cSWFhIZmZmYFlXq+X//qv/+K5555j3759OJ1OCgsLG2xXV1dHcXExTqcTAKfTSUFBQYMx9Z/rx5xqzpw5zJo1K/DZ5XKRkZHBxIkTcTgcjW7j8XjIyclhwoQJ2Gy25n/hEw5s3wjboIwYrvvWDS3eT2cQrDnfcPg9xhx/i6yarxg25cEgVtjxBGvOpWk036GnOQ+91pzz+iso5xLUAHP77bczfvz4BssmTZrE7bffzve//30AsrOzKS0tZcOGDYwY4e+s+sknn+Dz+Rg9enRgzC9+8Qs8Hk9gYnJycujXr1+jl48A7HY7drv9tOU2m+2ck9uUMWdTXeLvAVNiTSZe//I0yfnOefzFt8K/3mJg6Qq8nloio2ODWF3HdL5zLs2j+Q49zXnotcacN3V/zb6Jt6KigtzcXHJzcwHIy8sjNzeXAwcOkJyczIUXXtjgj81mw+l00q9fPwAGDBjANddcww9/+EPWrl3Lv//9b2bMmMEtt9wSeOT6tttuIyIigrvvvpstW7bw97//neeff77BGZa2pPpEF95ydeENmX4jx5NPF2It1Wxd/qbZ5YiISIg1O8CsX7+eYcOGMWzYMABmzZrFsGHDmDt3bpP3sWjRIvr378/VV1/NlClTuOyyyxr0eImPj2fZsmXk5eUxYsQI/uu//ou5c+e22UeoverCG3LWsDDyuk7yf9j8D3OLERGRkGv2JaRx48ZhGEaTx+/bt++0ZUlJSSxefPZHYC+66CJWrlzZ3PLMUd+FNzrN5EI6l9RLbod/vM6gii9wlR7HkaAzYCIinYXehRQEEdX+G4wtjsZvMJbW0XPQxey3ZmC3eNj+qXrCiIh0JgowQRB1oguvLf701yZI67FYrRzJ8Hd6jtquy0giIp2JAkwQxNcdByAm+QKTK+l8Mq/4TwAG1uRSdGS/ucWIiEjIKMCcJ8PnI8lXDEBclwyTq+l8Lug5gO3hAwizGOz+9M9mlyMiIiGiAHOeyl0lRJ7owpuUpvcgmaGsz40ApOx929xCREQkZBRgzlNp4UEAyo0oomLiTK6mc+p71R14jDB6e/ewf9sGs8sREZEQUIA5TxXH/T1gSq2NdwiW1pfYpStbYi4G4MjnfzG5GhERCQUFmPNUXXIEgPLwJJMr6dy8F/4HAD0Ov4/P6zW5GhERaW0KMOfJU+ZvYldjVxM1Mw0a910qjCi6cowd6z8yuxwREWllCjDnySj3N7HzRHUxuZLOLTI6lq2JVwLgWrvI5GpERKS1KcCcp7CqYwD4YhRgzBY94hYA+h//iNqaKpOrERGR1qQAc57sJ7rwhsXpNQJmG5A9lUKSiKeSrSv+aXY5IiLSihRgzlOM29+F156oAGO2sPBw9jonA2B8+XeTqxERkdakAHOeHN4SAKKT9BqBtqDLpXcAcGHFKspKikyuRkREWosCzHnweb0kGmUAOFL0Ise2oOegi8mzdifCUseOj9UTRkSko1KAOQ9lxYXYLP6eI4lddAamLbBYrRT0uB6AuB3/z+RqRESktSjAnIfSY4cAKCGOCHukydVIvV7j78ZrWBjg2cLhvVvMLkdERFqBAsx5qCjyd+Et02sE2pQu6T3YEjUcgAOfLTS3GBERaRUKMOehttQfYCpseo1AW+Me+B0AMg++g+HzmVyNiIgEmwLMeahz+bvw1thTTK5ETjXwqlupNCK5wChg+7ocs8sREZEgU4A5HxWFANTpNQJtTnRsPFsSrwKgfI2eRhIR6WgUYM5DeJU/wBCbam4h0qiYUd8DYMDxj6mpqjC5GhERCSYFmPMQeI2AQ11426IBY67hKF2Is1Sz5dO/mV2OiIgEkQLMeYj1+LvwRiYowLRF1rAw9nW7DgDb1wowIiIdiQLMeYj3FQMQk6wmdm1Vtyv+E4BB1espyj9gbjEiIhI0CjAt5HHXkmCUAxCvLrxtVkafIewI70+YxWD3x6+ZXY6IiASJAkwLlRYdxWoxqDOsJCTrElJbVtr3ZgDS9vxTPWFERDoIBZgWKjt2GIASSzxh4eEmVyNn03/893Eb4WT59rHn61VmlyMiIkGgANNClcX+LryuML1GoK2LT+rC147LATi+8k8mVyMiIsGgANNCtaVHAai0JZtciTSFbcTtAPQvWkpNdaXJ1YiIyPlSgGkhrysfgNpIvUagPRh02fUUkEw8leoJIyLSASjAtJC18sRrBKL1GoH2ICw8nL0XnOgJ89VfTa5GRETOlwJMC4VX+7vwWvQagXYj48ofAv6eMIWH80yuRkREzocCTAtFnniNQHi8HqFuL7r1vpCttgsJsxjs+ej3ZpcjIiLnQQGmheLq/F14oxLSTa5EmqNi4HcByNinnjAiIu2ZAkwLJfj870GKTVEX3vZk0Pg7qDLsdDOOsn1djtnliIhICynAtEBNdSUOqgBI0GsE2pWYuAQ2J14NQPnqheYWIyIiLdbsALNixQquu+460tPTsVgsvP3224F1Ho+HBx98kMGDBxMTE0N6ejp33HEHR44cabCP4uJipk2bhsPhICEhgbvvvpuKiooGY7766isuv/xyIiMjycjI4KmnnmrZN2wFJYWHAHAb4TgS1AemvYkdcycAg0o+oaqizORqRESkJZodYCorKxkyZAgvvfTSaeuqqqrYuHEjjzzyCBs3buSf//wnO3bs4Fvf+laDcdOmTWPLli3k5OSwZMkSVqxYwT333BNY73K5mDhxIt27d2fDhg08/fTTPProo7z66qst+IrB5zrxGoFiSyIWq05itTcDLp7IIUtXYiw1bM75i9nliIhICzT7JT6TJ09m8uTJja6Lj48nJ6fhfQW//e1vufjiizlw4ACZmZls27aNpUuXsm7dOkaOHAnAiy++yJQpU1iwYAHp6eksWrQIt9vNn/70JyIiIhg0aBC5ubk888wzDYKOWarqXyMQnoieQWp/LFYrB3vcRLe8l3BsXQw3/tTskkREpJla/S2EZWVlWCwWEhISAFi9ejUJCQmB8AIwfvx4rFYra9as4cYbb2T16tWMHTuWiIiIwJhJkybx5JNPUlJSQmLi6e8fqq2tpba2NvDZ5XIB/staHo+n0drql59p/Zmc/BqB5m7b2bV0zoMt66of4PnDK/T3bGXP12vI7D/c1HpaU1uZ885C8x16mvPQa805b+o+WzXA1NTU8OCDD3LrrbficDgAyM/PJzW1YfO38PBwkpKSyM/PD4zJyspqMCYtLS2wrrEAM3/+fObNm3fa8mXLlhEdHX3WOk89a3Qu1n1bASiui+SDDz5o1rbi19w5bw2p4cPI9q5n93u/YfPeaWaX0+rawpx3Jprv0NOch15rzHlVVVWTxrVagPF4PHznO9/BMAxefvnl1jpMwJw5c5g1a1bgs8vlIiMjg4kTJwbCU2M15uTkMGHCBGw2W5OPteHAP6EKolKzGDtlynnX3pm0dM5bw+bYWli5njHuzwm78lXsUTGm1tNa2tKcdwaa79DTnIdea855/RWUc2mVAFMfXvbv388nn3zSIEA4nU4KCwsbjK+rq6O4uBin0xkYU1BQ0GBM/ef6Maey2+3Y7fbTlttstnNOblPGNDhWjb8Lb5jDqX9ZWqi5c94aLrriJvJX/hwnx1j/2d8Y+a0fm1pPa2sLc96ZaL5DT3Meeq0x503dX9AfoakPL7t27eKjjz4iObnhY8bZ2dmUlpayYcOGwLJPPvkEn8/H6NGjA2NWrFjR4DpYTk4O/fr1a/TyUahFuY8DEJGgW3jbs7DwcPK63wxA9ObXTa5GRESao9kBpqKigtzcXHJzcwHIy8sjNzeXAwcO4PF4+Pa3v8369etZtGgRXq+X/Px88vPzcbvdAAwYMIBrrrmGH/7wh6xdu5Z///vfzJgxg1tuuYX0dH9b/ttuu42IiAjuvvtutmzZwt///neef/75BpeIzPTNawS6mlyJnK+eE+7Ba1gY6P6aAztzzS5HRESaqNkBZv369QwbNoxhw4YBMGvWLIYNG8bcuXM5fPgw7777LocOHWLo0KF07do18GfVqlWBfSxatIj+/ftz9dVXM2XKFC677LIGPV7i4+NZtmwZeXl5jBgxgv/6r/9i7ty5beIRasPnI9FXCkCcuvC2e2ndevF1zBgAjnz8isnViIhIUzX7Hphx48ZhGMYZ159tXb2kpCQWL1581jEXXXQRK1eubG55ra6yooxYi/9xbb1GoIMY8Z+wcjX9Ct6ntqYKe+TZn1oTERHzqY1sM9W/RqDSiCQmLsHcYiQoLhx7E4UkkYiLzR+fPViLiEjboADTTBXH/V14S6wJ5hYiQRNui2BPxk0A2L/6P5OrERGRplCAaabqE68RKA9PMrkSCaasCT/Ga1i4sDaX/ds3ml2OiIicgwJMM3nK/N2CqyNSTK5EgsmZ2YevYi4BID/neZOrERGRc1GAaSZfub+hnidKAaajibhsBgAXFX1A2fGCc4wWEREzKcA0kzWpO5vtQyFtkNmlSJANHHMNe8KyiLK42fb+b80uR0REzqLV30bd0Yy++X7gfrPLkFZgsVo5fuFd9PryEXrsXUyd5xHCbRHn3lBEREJOZ2BETnLRNXdTjAMnRXyZs8jsckRE5AwUYEROEhkVw85u/wFA9KZXzzFaRETMogAjcoreU+7DbYQxwLOVXZtWmF2OiIg0QgFG5BQp6d35Kv5KAMo+fcHkakREpDEKMCKNiL/yPgAuKvuEoiP7Ta5GREROpQAj0og+w8ay3TaQCIuXXUsWmF2OiIicQgFG5AxqLv4pAEMP/01nYURE2hgFGJEzGHL1LWwPH0CUxc2ef8w1uxwRETmJAozIGVisVnxX/xKA4UXvcWj3ZpMrEhGRegowImcxMHsyX0aOwmbxkv+OzsKIiLQVCjAi5xA75TEARpZ/zJ6vVplcjYiIgAKMyDn1uugS1sddDUDFBzoLIyLSFijAiDRB1xsfx2OEMaRmHVtWfWB2OSIinZ4CjEgTXNBzEBu7XA9A+Cfz8Hm9JlckItK5KcCINFGvb8+jyrDTr247a//6mNnliIh0agowIk2U4sxk8+AHARix60V2rP/E5IpERDovBRiRZhh10/1sjL0Cm8VL3Ps/oqykyOySREQ6JQUYkWawWK30/sFrHLGkkW4UsueP38fw+cwuS0Sk01GAEWkmR0IyFde9iscIY3jFCtb+4xmzSxIR6XQUYERaoO/wcWzo8zMAhmz+NXu+/sLkikREOhcFGJEWuvjWR/gychSRFg8J//gOezevMbskEZFOQwFGpIWsYWF0/8Hr7A7rRTJlpPy/G/VkkohIiCjAiJyHhBQnXWYsY3v4ABxU0u29W9ny7/fNLktEpMNTgBE5T/GJKWTct5TN9qHEWGrotexOvvz0TbPLEhHp0BRgRIIgJi6B3jM/IDc6m0iLh0Gf3cPq1x6kzuM2uzQRkQ5JAUYkSCKjYhg08x3WxU8i3OIje/8r7HnyMg7u/trs0kREOhwFGJEgskXYGXnf31g//ElcRNOvbgfJ/3c1a954Wg3vRESCSAFGJMgsVisjv/Vjqu5eyWb7UKIttYze+iu2/XqsnlISEQkSBRiRVuLM6M3ABz7hi77/TY1hY6D7a/otuZGNT1/L/h25ZpcnItKuNTvArFixguuuu4709HQsFgtvv/12g/WGYTB37ly6du1KVFQU48ePZ9euXQ3GFBcXM23aNBwOBwkJCdx9991UVFQ0GPPVV19x+eWXExkZSUZGBk899VTzv52IyaxhYYy57RFK717N2oQpeA0LwytXcsHiK1n7/DQO7vrS7BJFRNqlZgeYyspKhgwZwksvvdTo+qeeeooXXniBV155hTVr1hATE8OkSZOoqakJjJk2bRpbtmwhJyeHJUuWsGLFCu65557AepfLxcSJE+nevTsbNmzg6aef5tFHH+XVV19twVcUMZ8zsw8Xz/wrB2/5iE3RlxBu8XFxyRIyFo3lyycn8PWKt3SPjIhIM4Q3d4PJkyczefLkRtcZhsFzzz3Hww8/zPXXXw/AX/7yF9LS0nj77be55ZZb2LZtG0uXLmXdunWMHDkSgBdffJEpU6awYMEC0tPTWbRoEW63mz/96U9EREQwaNAgcnNzeeaZZxoEHZH2pseAkfQY8C+2r1lGzfJnuKjyC4ZUr4VP1rLvswwK+k6j3/jvk5DiNLtUEZE2rdkB5mzy8vLIz89n/PjxgWXx8fGMHj2a1atXc8stt7B69WoSEhIC4QVg/PjxWK1W1qxZw4033sjq1asZO3YsERERgTGTJk3iySefpKSkhMTExNOOXVtbS21tbeCzy+UCwOPx4PF4Gq23fvmZ1kvwac79eg2/EoZfycG9WzmS8wIXHVtCD99Bemz/Ne5tT7MpZgy+wd+l/+U3EWGPPK9jac5DS/Mdeprz0GvNOW/qPoMaYPLz8wFIS0trsDwtLS2wLj8/n9TU1IZFhIeTlJTUYExWVtZp+6hf11iAmT9/PvPmzTtt+bJly4iOjj5r3Tk5OWddL8GnOT9JxrdYmjoe68HPGVj+Of0s+xhW9W9Y829Kv/g5a23DOJY4knDnIMLCI869vzPQnIeW5jv0NOeh1xpzXlVV1aRxQQ0wZpozZw6zZs0KfHa5XGRkZDBx4kQcDkej23g8HnJycpgwYQI2my1UpXZqmvOz+TYAu7etp3DV/9E7/1+kWoq5sm4lHFtJVaGdbbGj8fa+hh4XTyUx9YIm7VVzHlqa79DTnIdea855/RWUcwlqgHE6/dftCwoK6Nq1a2B5QUEBQ4cODYwpLCxssF1dXR3FxcWB7Z1OJwUFBQ3G1H+uH3Mqu92O3W4/bbnNZjvn5DZljASX5vzMel+UTe+LsvHW1bFlzYeU575Fj2Of4rQUMaJyBXy5Ar78ObvCelPkvIz4wZPpPWzcOS81ac5DS/Mdeprz0GuNOW/q/oIaYLKysnA6nXz88ceBwOJyuVizZg333nsvANnZ2ZSWlrJhwwZGjBgBwCeffILP52P06NGBMb/4xS/weDyBL5KTk0O/fv0avXwk0hGFhYcz6NKpcOlUDJ+PXV9+TtH6f9AlfyW9vXvo491Nn8O74fBCqv5lZ0fkICq7ZpMw6Cp6DbkcW8TpgV5EpKNodoCpqKhg9+7dgc95eXnk5uaSlJREZmYmM2fO5Fe/+hV9+vQhKyuLRx55hPT0dG644QYABgwYwDXXXMMPf/hDXnnlFTweDzNmzOCWW24hPT0dgNtuu4158+Zx99138+CDD7J582aef/55nn322eB8a5F2xmK10mfYWPoMGwtAUf4B8r54D8veT+jlWkuixcXg2o2wbyPse4mqJXZ22vvh6jICe9YYPLV6RFtEOpZmB5j169dz5ZVXBj7X33dy5513snDhQh544AEqKyu55557KC0t5bLLLmPp0qVERn5zenvRokXMmDGDq6++GqvVys0338wLL7wQWB8fH8+yZcuYPn06I0aMICUlhblz5+oRapETUpyZpNwwHZiOz+slb8dGCr/6CNvBf9OzchMJlgoGub+Cw1/B4dcYDuzf9gSFjkH40oeT2Ceb7gNHYY88+w3uIiJtVbMDzLhx4zAM44zrLRYLjz32GI899tgZxyQlJbF48eKzHueiiy5i5cqVzS1PpNOxhoWRNXAUWQNHAeDzetm/M5eCLcuxHFqLs+xLMowjdDcO0b3sEJR9CNvA/U4Yu8N7UOwYgNF1CPE9R9B9wMVExcSZ/I1ERM6twzyFJCJ+1rAwug8YQfcB/nvMPB4P/+8ff6dXShg1+zcQfSyXzJrtJFrK6e3dAyV7oGQJbAXfexYOWrtyLKYPtSkDicoYSmqvYTgzemMNCzP5m4mIfEMBRqQTsEXFceEVU7DZbgXA8Pk4cmAXBTu+oObAJqKPb6ZbzU6SLWVkGEfIqDgCFcthH7ASKo1IDtm6UxbXG19Kf6LTB5Laawhp3XphseqdsCISegowIp2QxWolvUc/0nv0A+4MLC/KP8DRHeupPJBL+LEtJFfs5gLvQWIsNfSr2wElO6DkfdgFLIcqw85hWyZl0d3xJPYmIq0/iZkDSe85iMjoWNO+n4h0fAowIhKQ4swkxZkJ3BRY5nHXsn/vFor2fonnyNdElO4iqSqPC7xHiLbU0qduF7h2gesj2A+sBZ9hId+SwjF7NyrjemAk9iLK2YfEbv1xdu+rm4dF5LwpwIjIWdki7HTvP5zu/Yc3WO5x17I/bxvF+76kJn8H4cV7cFTm0bXuIA5LFU6O4aw9BrWboAj/WRvAa1g4aknhuP0CqmIy8CZkEZHSk/gL+tAlox/xSV1C/yVFpN1RgBGRFrFF2Onebyjd+w1tsNzw+SguOkph3hZcR3bgPbYLuyuP+OpDOOuOEGOpoSvH6Fp7DGpzoRjY+832LqI5FubEFZlObWwGlsTu2FN6EN+1J10y+hLrUDNLEVGAEZEgs1itJKVeQFLqBcDEBusMn4+iwsMU7d9K+dFd1B3Pw1a2n9iqQ3SpO0IyZTiowuHdC5V7oRJo+FYRyoihKCwNl91JbcwFkJBBRFImsWlZJKX3JKnLBXpiSqQTUIARkZCxWK2kODNIcWYAk05bX1VRRuHBXZQe2U1N4V4o3U9ExWEcNUdI8RaQQAXxVBLv3QtVe6EKONZwH24jnGPWZEptqVRHOvHEpmNN6IY9KYPYLpkkp/ckITlNT0+JtHMKMCLSZkTHxtNjwEgYMLLR9eVlxRQd2kVZfh61RfswSg5gqzhMTE0BiXWFpBjFRFjquMAo4AJ3Abi/BhdwpOF+ag0bRdYkysK7UB2ZiiemKzi6YkvoRkzyBTjSupPszNDNxiJtmAKMiLQbcfFJxMWPhkGjG13vcddScGQfpfl5VB3bj6fkEFbXYSKqjhJbW0iit4gUSrFbPP6Q4ykAD1AO5J++vxIclFiTqIjoQk1kF7yxTqxxTmwJ6UQnp5OQmklSWsY53wQuIsGnACMiHYYtwn5Sf5vGuWtrKDqyj7LC/VQVHcRTehhcR7BV5hNde4z4uiJSfMeJsNSRiItEnwtq9kENUNr4PkuJpdSaRLktmRp7ChHuCNa5NhKRkE5kYlfiki8gIbUbjsQuunQlEiQKMCLSqUTYI0nP6k96Vv8zjjF8PkqPF1BSeIDyYwepKTmCt+wI1vKjRFQfI9pdRHxdMUlGMREWLwlUkOCrgNoDUHtiJ3tO36/bCKPEkoArPIkqWxK19mS80V2wxKYSHp9GZIKTuKR0HCnpxCel6mZkkbNQgBEROYXFaiWhS1cSunQFGr9cBf4XZ5YUF1J67BAVRYeoKTmKp/QIriO7SLG5iXIfJ7aumERfMQ4qibB4SeM4aXXHoQ6o5oxndTxGGMctDlxhCVSFJ1JrT6IuKgUjOoXwuFRs8WlExacSl5xOfIqT6BiHzu5Ip6IAIyLSQtawMBK7dCWxS1fA/zZwj8fDBx98wJApU7DZbIGxtTVVlBQexlV0hKriI7jL8vGWF2CtPIat+hiR7mJi6kpI8JUQTyU2i5culNDFWwLePP+ZHdeZa6kxbJRaEigPi6falkitPRFvZDJGdDJhsV2IcKQSmZBKTGIajuSuOOKTFHikXVOAEREJAXtkNM7MPjgz+5xzrLu2htKio/6wU5JPbVkB3vJCqCgkrKaYiNrjRHtKiKsrIcEoI9LiIdLi8Xc/9h4DL/57dsrOfAyPEUaZJY5yazxV4fHURiTgsSfhi0rCEp1MeGwKEY4uRCWkEpOQSnxyms7ySJuiACMi0sZE2CNJvSCL1AuyzjnW8PmorHRRVnSU8uNHqS4rxO0qxFd+DCqPnQg8xUTVlRHrLSXe5yLGUoPN4iWFUlJ8peDG/6fi7MeqNWyUWeKotDpOhJ54PPZEfJGJWKKTCItJxhaXQmRcMtEJXYhLTMOR2IWwcP2qkeDTT5WISDtmsVqJiUsgJi4BsgY0aZuaqgrKiguoKC6gqvQYteWF1LmOYVQVY60+TnhtCXZ3CVF1ZcR5y4g3yrFbPNgtHlIpBl9xk0OPz7BQZomm3BJHVZiDmvB43BHxeO0JgeATHpOELS75m+CT0IW4hBQFHzkr/XSIiHQykdGxREbHktatV5PGGz4fVVXllB3Pp7L0GNVlx3C7iqirOIZRWYylpoTw2hJs7jKi6sqI8bqI87mIs1RjtRj+7slGJdTl+29ermlanS5iKLfEUmWNoybcgTsinrqIeHyRCViiEgmLTiQsJonw6ARqig9QeHgvCclOomPidKmrE1CAERGRs7JYrUTHxhMdGw/dz9xj51Qedy1lxYVUlh6jquwYta7juCuK8FUex6guxVpTiq22hAhPGVF1LqK95cQZFcRaqgFwUInDqARvgf++ntozH+sigIUPA/7H1cstsVRZYqgKi6M23IHb5sBrj8dnj8cSlUBYdCLhMQnYYpOIiksmypFMbEIKsXEJeny9nVCAERGRVmGLsJ/07qumCwSfsuNUu4qoLS/CXV6Mt7L4RPApIay2FJu7jMi6cqK85cT4ynGceHorwuIlmTKSjTL/GZ/6sz7l5z6298Qlr0pLLFXWWGrC43CHx1Fni8Nnj8eIjMcalUBYdAK22ETsMUnY4xKJdiQRG5+sG51DSAFGRETalOYGn/pH1ydfcw2V7mrKS4uoKjtOtes47opiPJUl+KpKoLoUS00pYe4ybB4X9rqKQPiJNSqJtHgIO/mSVxPO/JyqzrBSbomhyhJN9WkByIFhd2CJjMcaFY8tJhFbTDyRsYlExSUT40gkNj6JcFtEyyauk1GAERGRDqHBDc0ZvZu9fU11JRVlx6ksLaLadZzaylI8FcV4q0rxVZdgqS7F6i4n3OPC5in3n/3xVRJjVBJr+M/+hFt8JFJOolHeogAEUGXYqbDEUG2NocYaQ014HB5bLF5bHL6IODgRgsKiHIRHxxMRk4g9NoGouASi45KIjU/CFmFv9vdvbxRgREREgMioGCKjYkhxZjZ7W8Pno7qqgvLSIqrLi6kuL6G2ogRPZQne6jJ81aVQU47VXUa424XNU4G9rpxIXyXRvgpijCqiLf6kE22pJZpa/9NePpp143O9aiOCSks01ZZoqsNicYdF4w6Pw2uLxRsRBxFxEOnAGukgLDoeW1Q8ETHxRMYlEhmbQHRcIjGx8W36SbC2W5mIiEg7YbFaiYp1EBXrAHq2aB8edy0VZcVUnQhANRUleCrLqKsqxVtVilFbjqWm7MRZoHLC6yqw11WcCEGVDUJQlMVNFG4wSr+5D6iZZ4IAKo1IfxCyRlNrjaY2LBZPeCx1thgiq3xsjfcyZNxNLfq+50sBRkREpA2wRdhPejVFy9R53FS6SqgsLw2cBXJXlVFXWeo/E1TjgloXVncFYe5ywuoqiagrx+6tIvJECIox/O/tAoix1BBDTcOzQScFoVX7egMKMCIiInIewm0RxCenEZ+cdl77qamupKq8lCpXCTWVpdRUlFBXVU5ddRneahfe6jJKCg6Q2fuSIFXefAowIiIi0kD9/UBJqRc0ur7+ya/+oyeFuLJv6GF1ERERaXcUYERERKTdUYARERGRdkcBRkRERNodBRgRERFpdxRgREREpN1RgBEREZF2RwFGRERE2h0FGBEREWl3FGBERESk3Ql6gPF6vTzyyCNkZWURFRVFr169ePzxxzEMIzDGMAzmzp1L165diYqKYvz48ezatavBfoqLi5k2bRoOh4OEhATuvvtuKioqgl2uiIiItENBDzBPPvkkL7/8Mr/97W/Ztm0bTz75JE899RQvvvhiYMxTTz3FCy+8wCuvvMKaNWuIiYlh0qRJ1NTUBMZMmzaNLVu2kJOTw5IlS1ixYgX33HNPsMsVERGRdijoL3NctWoV119/PVOnTgWgR48e/PWvf2Xt2rWA/+zLc889x8MPP8z1118PwF/+8hfS0tJ4++23ueWWW9i2bRtLly5l3bp1jBw5EoAXX3yRKVOmsGDBAtLT04NdtoiIiLQjQQ8wl1xyCa+++io7d+6kb9++fPnll3z++ec888wzAOTl5ZGfn8/48eMD28THxzN69GhWr17NLbfcwurVq0lISAiEF4Dx48djtVpZs2YNN95442nHra2tpba2NvC5rKwM8F+K8ng8jdbq8Xioqqri+PHj2Gy2oHx/OTvNeehpzkNL8x16mvPQa805Ly8vB2hw60ljgh5gHnroIVwuF/379ycsLAyv18sTTzzBtGnTAMjPzwcgLS2twXZpaWmBdfn5+aSmpjYsNDycpKSkwJhTzZ8/n3nz5p22PCsr67y/k4iIiIRWeXk58fHxZ1wf9ADzxhtvsGjRIhYvXsygQYPIzc1l5syZpKenc+eddwb7cAFz5sxh1qxZgc8+n4/i4mKSk5OxWCyNbuNyucjIyODgwYM4HI5Wq02+oTkPPc15aGm+Q09zHnqtOeeGYVBeXn7O20WCHmBmz57NQw89xC233ALA4MGD2b9/P/Pnz+fOO+/E6XQCUFBQQNeuXQPbFRQUMHToUACcTieFhYUN9ltXV0dxcXFg+1PZ7XbsdnuDZQkJCU2q2eFw6Ic+xDTnoac5Dy3Nd+hpzkOvteb8bGde6gX9KaSqqiqs1oa7DQsLw+fzAf5LOk6nk48//jiw3uVysWbNGrKzswHIzs6mtLSUDRs2BMZ88skn+Hw+Ro8eHeySRUREpJ0J+hmY6667jieeeILMzEwGDRrEpk2beOaZZ7jrrrsAsFgszJw5k1/96lf06dOHrKwsHnnkEdLT07nhhhsAGDBgANdccw0//OEPeeWVV/B4PMyYMYNbbrlFTyCJiIhI8APMiy++yCOPPMJPfvITCgsLSU9P50c/+hFz584NjHnggQeorKzknnvuobS0lMsuu4ylS5cSGRkZGLNo0SJmzJjB1VdfjdVq5eabb+aFF14Iaq12u51f/vKXp116ktajOQ89zXloab5DT3Meem1hzi3GuZ5TEhEREWlj9C4kERERaXcUYERERKTdUYARERGRdkcBRkRERNqdThtgXnrpJXr06EFkZCSjR48OvGxSgm/+/PmMGjWKuLg4UlNTueGGG9ixY4fZZXUqv/71rwMtDKT1HD58mO9973skJycTFRXF4MGDWb9+vdlldVher5dHHnmErKwsoqKi6NWrF48//vg536EjTbdixQquu+460tPTsVgsvP322w3WG4bB3Llz6dq1K1FRUYwfP55du3aFpLZOGWD+/ve/M2vWLH75y1+yceNGhgwZwqRJk07r/ivBsXz5cqZPn84XX3xBTk4OHo+HiRMnUllZaXZpncK6dev43//9Xy666CKzS+nQSkpKuPTSS7HZbPzrX/9i69at/OY3vyExMdHs0jqsJ598kpdffpnf/va3bNu2jSeffJKnnnqKF1980ezSOozKykqGDBnCSy+91Oj6p556ihdeeIFXXnmFNWvWEBMTw6RJk6ipqWn94oxO6OKLLzamT58e+Oz1eo309HRj/vz5JlbVeRQWFhqAsXz5crNL6fDKy8uNPn36GDk5OcYVV1xh3HfffWaX1GE9+OCDxmWXXWZ2GZ3K1KlTjbvuuqvBsptuusmYNm2aSRV1bIDx1ltvBT77fD7D6XQaTz/9dGBZaWmpYbfbjb/+9a+tXk+nOwPjdrvZsGED48ePDyyzWq2MHz+e1atXm1hZ51FWVgZAUlKSyZV0fNOnT2fq1KkNft6ldbz77ruMHDmS//iP/yA1NZVhw4bx+9//3uyyOrRLLrmEjz/+mJ07dwLw5Zdf8vnnnzN58mSTK+sc8vLyyM/Pb/Dfl/j4eEaPHh2S36dB78Tb1hUVFeH1eklLS2uwPC0tje3bt5tUVefh8/mYOXMml156KRdeeKHZ5XRof/vb39i4cSPr1q0zu5ROYe/evbz88svMmjWLn//856xbt46f/exnREREcOedd5pdXof00EMP4XK56N+/P2FhYXi9Xp544gmmTZtmdmmdQn5+PkCjv0/r17WmThdgxFzTp09n8+bNfP7552aX0qEdPHiQ++67j5ycnAav6JDW4/P5GDlyJP/zP/8DwLBhw9i8eTOvvPKKAkwreeONN1i0aBGLFy9m0KBB5ObmMnPmTNLT0zXnnUCnu4SUkpJCWFgYBQUFDZYXFBTgdDpNqqpzmDFjBkuWLOHTTz+lW7duZpfToW3YsIHCwkKGDx9OeHg44eHhLF++nBdeeIHw8HC8Xq/ZJXY4Xbt2ZeDAgQ2WDRgwgAMHDphUUcc3e/ZsHnroIW655RYGDx7M7bffzv3338/8+fPNLq1TqP+dadbv004XYCIiIhgxYgQff/xxYJnP5+Pjjz8mOzvbxMo6LsMwmDFjBm+99RaffPIJWVlZZpfU4V199dV8/fXX5ObmBv6MHDmSadOmkZubS1hYmNkldjiXXnrpae0Bdu7cSffu3U2qqOOrqqrCam34aywsLAyfz2dSRZ1LVlYWTqezwe9Tl8vFmjVrQvL7tFNeQpo1axZ33nknI0eO5OKLL+a5556jsrKS73//+2aX1iFNnz6dxYsX88477xAXFxe4NhofH09UVJTJ1XVMcXFxp91jFBMTQ3Jysu49aiX3338/l1xyCf/zP//Dd77zHdauXcurr77Kq6++anZpHdZ1113HE088QWZmJoMGDWLTpk0888wz3HXXXWaX1mFUVFSwe/fuwOe8vDxyc3NJSkoiMzOTmTNn8qtf/Yo+ffqQlZXFI488Qnp6OjfccEPrF9fqzzm1US+++KKRmZlpREREGBdffLHxxRdfmF1ShwU0+ue1114zu7RORY9Rt7733nvPuPDCCw273W7079/fePXVV80uqUNzuVzGfffdZ2RmZhqRkZFGz549jV/84hdGbW2t2aV1GJ9++mmj//2+8847DcPwP0r9yCOPGGlpaYbdbjeuvvpqY8eOHSGpzWIYalkoIiIi7UunuwdGRERE2j8FGBEREWl3FGBERESk3VGAERERkXZHAUZERETaHQUYERERaXcUYERERKTdUYARERGRdkcBRkRERNodBRgRERFpdxRgREREpN1RgBEREZF25/8Ds2d98BY66IoAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "nb_pts = 100\n", "\n", "temp = 400\n", "pres = 200000\n", "phi = np.linspace(0.3,10,nb_pts)\n", "\n", "Tad_gas1 = np.zeros(nb_pts)\n", "Tad_gas2 = np.zeros(nb_pts)\n", "\n", "gas1 = ct.Solution('Mechanisms/GRI-Mech2.11.yaml')\n", "gas2 = ct.Solution('gri30.yaml')\n", "\n", "for i in range(0,len(phi)):\n", " gas1.TP = temp,pres\n", " gas2.TP = temp,pres\n", " gas1.set_equivalence_ratio(phi[i],'CH4:1',{'O2':1,'N2':3.76})\n", " gas2.set_equivalence_ratio(phi[i],'CH4:1',{'O2':1,'N2':3.76})\n", " gas1.equilibrate('HP')\n", " gas2.equilibrate('HP')\n", " Tad_gas1[i] = gas1.T\n", " Tad_gas2[i] = gas2.T\n", "plt.figure\n", "plt.plot(phi,Tad_gas1,label='GRI-Mech2.11')\n", "plt.plot(phi,Tad_gas2,label='gri30')\n", "plt.grid()\n", "plt.legend()\n", "\n", "# L'image présentée n'est pas pour de méthane mais de l'éthylène (C2H4)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 4. Compute the low heating value and high heating value of a fuel" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In combustion, the low heating value (LHV or PCI in French) corresponds to the energy released when the fuel and $O_2$ are converted into $CO_2$ and $H_2O$. It is directly linked to the formation enthalpy of the species involved and given by :\n", "\n", "$$ LHV = \\sum_i^c n_i \\bar{h_i}(T_ref) $$\n", "\n", "Here, we want to compute the LHV for methane. The complete reaction can be written as:\n", "\n", "$$ CH_4 + 2 O_2 \\to CO_2 + 2 H_2O $$\n", "\n", "To compute the LHV, we will use the thermodynamics data from the gri30 mechanism. The mixture will be fixed at the reference state (298K, 1 atm). \n", "\n", "First, we need to compute the enthalpy of the reactants : " ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "import cantera as ct" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", " gri30:\n", "\n", " temperature 298 K\n", " pressure 1.0132e+05 Pa\n", " density 1.0911 kg/m^3\n", " mean mol. weight 26.68 kg/kmol\n", " phase of matter gas\n", "\n", " 1 kg 1 kmol \n", " --------------- ---------------\n", " enthalpy -9.3222e+05 -2.4871e+07 J\n", " internal energy -1.0251e+06 -2.7349e+07 J\n", " entropy 7652.5 2.0417e+05 J/K\n", " Gibbs function -3.2127e+06 -8.5713e+07 J\n", " heat capacity c_p 1179.9 31480 J/K\n", " heat capacity c_v 868.29 23166 J/K\n", "\n", " mass frac. Y mole frac. X chem. pot. / RT\n", " --------------- --------------- ---------------\n", " O2 0.79956 0.66667 -25.079\n", " CH4 0.20044 0.33333 -53.622\n", " [ +51 minor] 0 0 \n", "\n" ] } ], "source": [ "gas = ct.Solution('gri30.yaml')\n", "\n", "# Set the reactants state\n", "gas.TP = 298, 101325\n", "gas.set_equivalence_ratio(1.0, 'CH4:1.0', 'O2:1.0')\n", "\n", "gas()\n", "\n", "# Get the enthalpy and the mass fractions of the fuel\n", "h_reactants = gas.enthalpy_mass\n", "Y_CH4 = gas[\"CH4\"].Y[0]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now that we have the enthalpy of the reactants, we need to compute the enthalpy of the products :" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", " gri30:\n", "\n", " temperature 298 K\n", " pressure 1.0132e+05 Pa\n", " density 1.0911 kg/m^3\n", " mean mol. weight 26.68 kg/kmol\n", " phase of matter gas\n", "\n", " 1 kg 1 kmol \n", " --------------- ---------------\n", " enthalpy -1.0959e+07 -2.9239e+08 J\n", " internal energy -1.1052e+07 -2.9487e+08 J\n", " entropy 7587.1 2.0242e+05 J/K\n", " Gibbs function -1.322e+07 -3.5271e+08 J\n", " heat capacity c_p 1303.1 34767 J/K\n", " heat capacity c_v 991.5 26453 J/K\n", "\n", " mass frac. Y mole frac. X chem. pot. / RT\n", " --------------- --------------- ---------------\n", " H2O 0.45016 0.66667 -120.72\n", " CO2 0.54984 0.33333 -185.63\n", " [ +51 minor] 0 0 \n", "\n" ] } ], "source": [ "# Set the products state\n", "X_products = {'CO2': gas.elemental_mole_fraction('C'), 'H2O': 0.5*gas.elemental_mole_fraction('H')}\n", "gas.TPX = 298, 101325, X_products\n", "#or\n", "#gas.equilibrate('TP') (sometimes you will get minor species that can alter a bit the LHV calculation)\n", "\n", "gas()\n", "\n", "# Get the products enthalpy\n", "h_products = gas.enthalpy_mass" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, we can compute the LHV for methane :" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "LHV [MJ/kg] = 50.02548811553696\n" ] } ], "source": [ "LHV = -(h_products - h_reactants) /Y_CH4\n", "print(\"LHV [MJ/kg] = \", LHV/1e6)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Conclusion" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We have just generated the skeleton of a script to perform a series of common equilibrium calculations\n", "to obtain the constant pressure equilibrium composition of a fuel/air mixture. Starting from there,\n", "you could modify your initial conditions, plot the mole/mass fractions of other species, change the\n", "solver or even try another fuel (methane, acetylene) without changing your mechanism.
\n", "Technically, adiabatic flame calculations could also be performed at constant volume: simply invoke\n", "the good equilibrate option of your equilibrate function, 'UV' (see 3.1.2), in your script.\n", "The computation of the low heating value for a fuel has also been investigated in this tutorial" ] } ], "metadata": { "kernelspec": { "display_name": "cantera-avbp3.0-py39", "language": "python", "name": "cantera-avbp3.0-py39" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.5" }, "vscode": { "interpreter": { "hash": "e6a4b3b1c4ad517cd8271eeabb831a2b1c3be73abb21c806c4703ebe5f8b2686" } } }, "nbformat": 4, "nbformat_minor": 2 }