09Air_flow_ventilation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "385f59fc-dfa8-4336-9b92-8990fa583f7e",
   "metadata": {},
   "source": [
    "# [Weather data and solar radiation](01WeatherData.ipynb)\n",
    "\n",
    "# [Thermal circuit and state-space representation for a thermal circuit with capacities in every node: simple wall](02SimpleWall.ipynb)\n",
    "\n",
    "# [Thermal circuit and state-space representation for a thermal circuit with capacities in some nodes: cubic building](03CubicBuilding.ipynb)\n",
    "\n",
    "# [Thermal circuits assembling](04AssemblingCircuits.ipynb)\n",
    "\n",
    "# [Switch between models: heating & cooling and free-running](05SwitchModels.ipynb)\n",
    "\n",
    "# [Control input: heating & cooling and free-running](06Control_Input.ipynb)\n",
    "\n",
    "# [Radiation coupled with convection](07Coupled_rad_convection.ipynb)\n",
    "\n",
    "# [Sensible thermal load in steady-state](08Thermal_load.ipynb)\n",
    "\n",
    "# Air flow by ventilation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5aaefbfc-f6c2-49b5-83d0-06eac3fa1635",
   "metadata": {},
   "source": [
    "## Problem\n",
    "The natural ventilation system of a two-storey building is composed of (Figure 1):\n",
    "- Orrifices located on the north and south faces located at hight $z_0 = 0.5 \\, \\mathrm{m}$. Every façade has two orrifices with the pressure exponent $n_o = 0.64$ and the flow coefficient $K_o = 10 \\, \\mathrm{m^3/(h\\,Pa}^{n_o})$.\n",
    "- An exhaust air chimney on the south side of the roof at hight $z_c = 8.5 \\, \\mathrm{m}$. This chimney may be considered an orrifice with the pressure exponent $n_c = 0.50$ and the flow coefficient $K_c = 50 \\, \\mathrm{m^3/(h\\,Pa}^{n_c})$.\n",
    "- The other elements of the building are considered airthight.\n",
    "\n",
    "![vent_bldg](./figures/09_nat_vent.svg)\n",
    "> Figure 1. Natural ventilation with wind and stack effect.\n",
    "\n",
    "The static pressure at the level of the orrifices is considered $P_o = 0 \\, \\mathrm{Pa}$.\n",
    "\n",
    "Let's consider that the outdoor air has the parameters $\\theta_o = -5 \\, \\mathrm{°C}$, $\\varphi_o = 100 \\, \\%$,$\\rho_o = 1.315 \\, \\mathrm{kg/m^3}$ and the indoor air has the parameters $\\theta_i = 22 \\, \\mathrm{°C}$, $\\varphi_i = 50 \\, \\%$, $\\rho_i = 1.191 \\, \\mathrm{kg/m^3}.$\n",
    "\n",
    "In these conditions, estimate:\n",
    "- The possible domain of variation for the indoor pressure $p_i$ for the flow configuration shown in Figure 1.\n",
    "- The value of the indoor pressure $p_i$.\n",
    "- The mass flow rates:\n",
    "    - $\\dot m_{ni}$: from outdoors to indoors through the nord façade;\n",
    "    - $\\dot m_{ci}$: from outdoors to indoors through the south façade;\n",
    "    - $\\dot m_{ic}$: from indoors to outdoors through the chimney;\n",
    "    \n",
    "in two situations:\n",
    "1. Neglecting the wind effect.\n",
    "2. Considering that the wind blows from North to South with speed $v = 18 \\, \\mathrm{km/h}$ and that the pressure coefficients are $C_{p,n} = 0.4$ for the north face, $C_{p,s} = -0.55$ for the south face, and $C_{p,c} = -0.7$ for the chimney. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "45299720-7f28-40e8-b661-44a2f50f9eb6",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from scipy.optimize import fsolve\n",
    "\n",
    "# Data\n",
    "# standard acceleration due to gravity, m/s²\n",
    "g = 9.8\n",
    "# pressure exponent, dimensionless\n",
    "no, nc = 0.64, 0.50             # orrifice, cheminée\n",
    "# flow coefficients, m³/(s·Paⁿ)\n",
    "Ko, Kc = 10 / 3600, 50 / 3600   # orrifice, chimney\n",
    "# density, kg/m³\n",
    "ρi, ρo = 1.191, 1.315           # indoor, outdoor\n",
    "# pressure coefficients, dimensioneless\n",
    "Cpn, Cps, Cpc = 0.4, -0.55, -0.7    # north, south, chimney\n",
    "# hignt, m\n",
    "z = 8"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "83d8de71-fdce-4947-b216-6b8376d8eeed",
   "metadata": {},
   "source": [
    "## Summary of theory\n",
    "\n",
    "### Sources\n",
    "#### Pressure sources\n",
    "\n",
    "Pressure sources (Rock et al. 2017):\n",
    "- stack pressure: $P_s = P_0 - \\rho g z$;\n",
    "- wind pressure: $P_w = C_p \\frac{1}{2} \\rho_o v^2$\n",
    "\n",
    "where:\n",
    "- $P_s$ stack pressure, Pa\n",
    "- $P_0$ stack pressure at reference height, Pa\n",
    "- $g$ gravitational acceleration, m/s²\n",
    "- $\\rho$ indoor or outdoor air density, kg/m³\n",
    "- $z$ height above the reference level, m\n",
    "- $P_w$ wind pressure difference relative to the outdoor pressure in undisturbed flow, Pa\n",
    "- $C_p$ pressure coefficient, dimensionless\n",
    "- $\\rho_o$ outside air density, kg/m³\n",
    "- $v$ wind speed, m/s\n",
    "\n",
    "#### Flow sources\n",
    "Flow sources inject a controlled volumetric airflow rate into the control volume, as is the case in [HVAC systems](https://en.m.wikipedia.org/wiki/Ventilation_(architecture)).\n",
    "\n",
    "### Constitutive equation\n",
    "\n",
    "Each opening of a building can be modelled by a [constitutive equation](https://en.m.wikipedia.org/wiki/Constitutive_equation) or power law (Rock et al. 2017):\n",
    "\n",
    "$$ \\dot V_{i,j} = K_{i,j} \\Delta P^n_{i,j} = K_{i,j}(P_i - P_j)^n$$\n",
    "\n",
    "where:\n",
    "- $\\dot V_{i,j}$ volumetric airflow through opening from control volume $i$ to control volume $j$, m³/s\n",
    "- $n$ pressure exponent, dimensionless; typical value about 0.65\n",
    "- $K_{i,j}$ flow coefficient, m³/(s·Paⁿ)\n",
    "\n",
    "The values of the pressure exponent $n$ and flow coefficient $K_{i,j}$ cann be obtained experimentally by using the [blower door](https://en.wikipedia.org/wiki/Blower_door) test.\n",
    "\n",
    "### Mass conservation equation\n",
    "\n",
    "The [mass conservation](https://en.m.wikipedia.org/wiki/Conservation_of_mass) law in steady-state states that the algebraic sum of mass flow rates entering in a control volume $j$ is zero:\n",
    "\n",
    "$$\\sum_i \\dot m_{i,j} = 0$$\n",
    "\n",
    "Mass air flow rate from control volume $i$ to control volume $j$ is:\n",
    "\n",
    "$$ \\dot m_{i,j} = \\rho_i \\, \\dot V_{i,j}$$\n",
    "\n",
    "where:\n",
    "- $\\dot m_{i,j}$ mass flow rate of air from control volume $i$ to control volume $j$, kg/s \n",
    "- $\\rho_i$ air density in control volume $i$, kg/m³\n",
    "- $\\dot V_{i,j}$ volumetric air flow rate from control volume $i$ to control volume $j$, m³/s\n",
    "\n",
    "**Notes:**\n",
    "1. The constitutive law relates volumetric airflow to pressure difference.\n",
    "2. The balance equation is done for mass flow rate.\n",
    "3. Generally, $\\dot m_{i,j} \\neq \\dot m_{j,i}$, even if $\\dot V_{i,j} = \\dot V_{j,i}$. This implies that if the mass flow rates obtained by slvinf the aeraulic circuit are negative, than the direction of the flows need to changed and a new solution to be searched for.\n",
    "\n",
    "### Domain of variation for the indoor pressure\n",
    "\n",
    "For the flow directions shown in Figure 1,\n",
    "\n",
    "$$\\dot m_{i,c} > 0$$\n",
    "\n",
    "the indoor pressure $p_i$ needs to satisfy the inequality (see Figure 2a):\n",
    "\n",
    "$$p_i + [P_{si} - (P_{dc} + P{sc})] > 0$$\n",
    "\n",
    "or the inequality (see Figure 2b):\n",
    "\n",
    "$$p_i > P_{dc} + P{sc} - P_{si}$$\n",
    "\n",
    "![vent_bldg](./figures/09_nat_vent_1P.svg)\n",
    "> Figure 2. The model from Figure 1 with only one pressure source on a branch: a) pressure source before conductance; b) pressure source after conductance."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "715aee41-50e8-47f8-8740-64921d122827",
   "metadata": {},
   "source": [
    "## Static pressure\n",
    "\n",
    "- $P_{si} = - \\rho_i g z$ indoors\n",
    "- $P_{so} = - \\rho_o g z$ outdoors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "e9e6af66-0b0a-4330-8493-82d738a871d1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Static pressure, Pa\n",
    "Psi = -ρi * g * z   # indoors\n",
    "Pso = -ρo * g * z   # outdoors"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9f383a27-1523-4636-a6d1-b024f807261c",
   "metadata": {},
   "source": [
    "## Neglecting the wind effect\n",
    "\n",
    "### Dynamic pressure\n",
    "\n",
    "The dynamic pressures are:\n",
    "- $P_{dn} = C_pn \\frac{1}{2} \\rho_o v^2$ - north face\n",
    "- $P_{ds} = C_ps \\frac{1}{2} \\rho_o v^2$ - south face\n",
    "- $P_{dc} = C_pc \\frac{1}{2} \\rho_o v^2$ - chimney"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "79547d54-dd96-498d-87d4-7039e32061c1",
   "metadata": {},
   "outputs": [],
   "source": [
    "v = 0       # m/s, wind speed\n",
    "# dynamic pressures, Pa\n",
    "Pdn = Cpn * 0.5 * ρo * v**2     # north face\n",
    "Pds = Cps * 0.5 * ρo * v**2     # south face\n",
    "Pdc = Cpc * 0.5 * ρo * v**2     # chimney"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38f777e0-04ed-4c3a-9e19-24416693b85d",
   "metadata": {},
   "source": [
    "### Domain of variation of indoor pressure\n",
    "The domain of variation of $p_i$ in order to have the flows shown in Figure 1 is\n",
    "\n",
    "$$P_{si} - (P_{dc} + P{sc}) \\le p_i \\le \\min (P_{dn},P_{ds})$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "bac36db5-e614-4b19-95d9-092e351a1ce2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Pressure domain: -9.72 Pa < pi < 0.00 Pa\n"
     ]
    }
   ],
   "source": [
    "pi_min, pi_max = -(Psi - Pdc - Pso), min(Pdn, Pds)\n",
    "print(f\"Pressure domain: {pi_min:.2f} Pa < pi < {pi_max:.2f} Pa\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a765e254-f5de-4ac5-850f-c4f435377d27",
   "metadata": {},
   "source": [
    "### Graphical solution of nonlinear equation for indoor pressure\n",
    "\n",
    "The mass flow balance equation for pressure node $p_i$ is:\n",
    "$$ \\dot m_{ni} + \\dot m_{si} - \\dot m_{ic} = 0$$\n",
    "\n",
    "where:\n",
    "- $ \\dot m_{ni} = 2 \\rho_o K_o (P_{dn} - P_{i})^{n_o}$ - mass flow through the two orifices on north façade;\n",
    "- $ \\dot m_{si} = 2 \\rho_o K_o (P_{ds} - P_{i})^{n_o}$ - mass flow through the two orifices on south façade;\n",
    "- $ \\dot m_{ic} = \\rho_i K_c (p_i + P_{si} - P_{dc} - P_{so})^{n_c}$ - mass flow from indoor air through the chimney.\n",
    "\n",
    "The mass balance is satisfied for the indoor pressure $p_i$ which solves the equation:\n",
    "\n",
    "$$f(p_i) = 0$$\n",
    "\n",
    "where:\n",
    "\n",
    "$$f(p_i) = 2 \\rho_o K_o (P_{dn} - P_{i})^{n_o} + 2 \\rho_o K_o (P_{ds} - P_{i})^{n_o} - \\rho_i K_c (p_i + P_{si} - P_{dc} - P_{so})^{n_c}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "ddd0b017-e982-4c8c-aea9-6a3afc72d086",
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(pi):\n",
    "    \"\"\"\n",
    "    Error in mass balance equation as a function of indoor pressure.\n",
    "    Used to find indoor pressure pi which makes the mass balance zero.\n",
    "    The north and south flows are enteriing.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    pi : float\n",
    "        Indoor pressure, Pa.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    y : float\n",
    "        Mass flow rate unbalanced in the mass balance equation\n",
    "    \"\"\"\n",
    "    y = 2 * ρo * Ko * (Pdn - pi)**no \\\n",
    "        + 2 * ρo * Ko * (Pds - pi)**no \\\n",
    "        - ρi * Kc * (pi + Psi - Pdc - Pso)**nc\n",
    "    return y"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ab182701-b75f-46a9-9cee-49ec1a61b390",
   "metadata": {},
   "source": [
    "The solution of this non-linear equation can be found graphically."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "0327f689-68e2-4bf0-b2b6-cc374e11327e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZYAAAEKCAYAAAAxXHOuAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAyw0lEQVR4nO3debyV0/7A8c9XM1I00UDRQKWORmU6CrcyRKG6Vxlu0oSQFCFDZZYxJVHhppuim26GdEyJBkejo1OGDu41/ByUQvn+/ljPubbTPufsvXuePX7fr9d6nf1Ma3+X4XzPs55nrSWqijHGGOOXfRIdgDHGmPRiicUYY4yvLLEYY4zxlSUWY4wxvrLEYowxxlflEx1AotWsWVMbNmwY9tj27dvZb7/9wh7Ly3M/mzULKLAEKq3d6S5T256p7YbMbXteHuzevZvmzcvFdP2qVau+VdVa4Y5lfGJp2LAhK1euDHssJyeH7OzssMeKdufkBBJWQpXW7nSXqW3P1HZD5rY9OxsKCwtZubJ6TNeLyGclHbOuMGOMMb7K+DuWWI0dm+gIjDEmdmPHwocffgZU971uSywxOuWUREdgjDGxO+UUKF/++0Dqtq6wGOXmumKMMakoNxfy8/cPpG67Y4nRiBHuZzo+vDfGpL8RI6CwsDEDB/pft92xGGOM8ZUlFmOMMb6yrrAY/fQTfPMNqIJIoqMxxpjkYXcsMdq+HbZutQf4xhhTnCWWGN17r7tTmTcv0ZEYY0z0JkyAgQO3BFJ30iUWEekmInkiki8io8McFxF50Du+RkTahByrLiJzReQjEdkoIp2CivP00+HEE2H+/KC+wRhjgtO5M7Rs+WMgdSdVYhGRcsAjQHegOdBPRJoXO6070MQrg4DJIcceABar6pFAa2BjULEuWwbHHAPr18PHHwf1LcYYE4xly2DdugMCqTupEgvQAchX1S2q+iswG+hZ7JyewEx1lgPVReQQETkAOBF4AkBVf1XVwqACvf56WL7cfba7FmNMqrn+epg27fBA6k62t8LqAVtDtguAjhGcUw/YBXwDPCkirYFVwJWqur34l4jIINzdDnXq1CGnhFGO27ZtK/FYYWEWAE2blmPGDKVjx9WlNiyVlNbudJepbc/UdkPmtr2wMIvdu3cH0vZkSyzhXtzVCM8pD7QBLlfV90TkAWA0cOMeJ6tOBaYCtGvXTkuaMru06bSrV3c/zz8fbrgBmjTJpl69sKemnEydRhwyt+2Z2m7I3LZXr+6mzQ+i7cnWFVYANAjZrg98GeE5BUCBqr7n7Z+LSzSBOucc9/OFF4L+JmOMSQ3JllhWAE1EpJGIVAT6AguKnbMAGOC9HXYs8IOqfqWq/wG2ikjRmo5dgQ1BB3zUUW4VSXvOYowxTlJ1hanqLhEZDrwMlAOmq+p6ERnsHX8MWAT0APKBn4GLQ6q4HHjGS0pbih3z1aRJf3zu1Qvuugu++IK06Q4zxqS3SZNg5cp8oJ3vdSdVYgFQ1UW45BG677GQzwoMK+HaXIL4pxRGVtYfny+5xP1LGjgQFi2yKV6MMckvKwsKC7cFUneydYWljNdecwWgcWO4+25YvBgee6z064wxJhm89hqsWnVgIHVbYonR7be7UmToUPjLX+Caa2zApDEm+d1+O8yadVggdVti8YkITJ8OVarABRfAb78lOiJjjEkMSyw+qlsXpkyBFSvcHYwWH4FjjDEZwBKLz8491w2YnDYN7rwz0dEYY0z8Jd1bYengtttgyxYYMwYaNYI+fRIdkTHGxI8llhhNmVLyMRF48km3ENiFF0Lt2nDyyfGLzRhjyjJlCrz3Xh57Tse496wrLEbNmrlSkkqV3DQvRxwBZ5wBb74Zt9CMMaZMzZrBoYfuCKRuSywx+te/XClNjRrw+utw6KHQowe8/XZ8YjPGmLL861+wbFmNQOq2xBKje+91pSx16rjkUq8edO8Ob70VfGzGGFOWe++FOXMalH1iDCyxxMEhh8DSpS65nHYaLFyY6IiMMSY4lljipG5dd7fSsiWcfTbMnJnoiIwxJhiWWOKoVi3XLZad7d4Wu+suG0RpjEk/Zb5uLCIHRVDP70GuL59OqlaFl15yieW662DTJnj0UahQIdGRGWOMPyIZx/KlV0qbDL4ccKgvEaWIWbNiv7ZSJXj2WTcr8vjx8Omn8M9//rHcsTHGBG3WLHj33Y1AJ9/rjiSxbFTVY0o7QUQ+8CmelNFgL1+m2GcfN7to48YwaBB07AgvvghHHulPfMYYU5oGDWDz5l8CqTuSZyyRpDP/U16Se+45V/bWRRfBkiVQWAgdOpQ9NsYYY/zw3HPw+uu1Aqm7zMSiqjsBROQ8Eanqfb5RROaJSJvQczLJ5Mmu+OGEE2DlSmjaFM46C265BX7/3Z+6jTEmnMmTYcGCYNZSj+atsBtV9ScROR44DZgB+PSr1TRo4F5HHjAAxo1zI/W//TbRURljTPSiSSy7vZ+nA5NV9UWgov8hZa4qVeCpp9zkcEuXQps2sHx5oqMyxpjoRJNYvhCRKcD5wCIRqRTl9SYCIu5h/rJlUK6c6ya76y7rGjPGpI4yE4OIdBIRwSWUl4Fu3piVg4Brgw0vc7VtCx98AD17uvEuPXrA118nOipjjClbJHccFwKrgOnAAcBPAKr6laq+EmBsSW3uXFeCVL26G9/y2GPwxhvQqhUsXhzsdxpjMsPcuXDLLesDqTuSt8IGq2obYBxwIPCUiLwrIhNE5EQRKRdIZEmuZk1XgiYCl10G77/vpoTp3h2uvBJ2BLOMgjEmQ9SsCdWq/RZI3RE/I1HVj1T1flXtBnQB3gbOA94LJLIk99RTrsTL0UfDihUuqTz4ILRr57rKjDEmFk89BYsXHxxI3TE9fFfVHaq6SFUvV9V2fgYkIt1EJE9E8kVkdJjjIiIPesfXFI2lCTleTkQ+EJFAJ6ePd2IBqFwZJk1y3WHff+8GVI4fD7t2xTcOY0zqCzKxRLzmvYhcHWb3D8AqVc31IxivW+0R4FSgAFghIgtUdUPIad2BJl7piBtLE7po85XARtzzoLT0l7/AunUwbBiMHQsLFrj/SI46KtGRGWNMdHcs7YDBQD2vDAKygcdFZJRP8XQA8lV1i6r+CswGehY7pycwU53lQHUROQRAROrjxtlM8ymepHXQQfCPf8Ds2ZCfD8ccA3ffDbt3l32tMcYEKeI7FqAG0EZVtwGIyM3AXOBE3Ftjd/kQTz1ga8h2AX++GynpnHrAV8AkYBRQtbQvEZFBuMRInTp1yMnJCXvetm3bSjxWWJgFQE5ObmlfFbg6deDxxyswaVJTRo2qxVNP/cCoUXkcdtjPMddZWrvTXaa2PVPbDZnb9sLCLHbv3h1I26NJLIcCv4Zs/wYcpqo7RMSvKTLDTc1ffCmssOeIyBnA16q6SkSyS/sSVZ0KTAVo166dZmeHPz0nJ4eSjhVNcV/S8Xg75xx393L55dUYNKgDN90Eo0bFts5Lae1Od5na9kxtN2Ru26tXh8LCwkDaHk1X2LPAchG52btbeQf4h4jsB2wo/dKIFQChE9LXx60FE8k5xwFnicinuC60LiLytE9x7WHRIleShQj06wcbNrgkM3YstG/vJrc0xpjiFi2CO+5YE0jd0SSWRcClQCHuof1g4ANV3a6qf/MpnhVAExFpJCIVgb7AgmLnLAAGeG+HHQv84A3WHKOq9VW1oXfd66p6gU9x7WHffV1JNrVruzuX+fPdSP2OHeGaa2D79kRHZoxJJvvuC5UrBzNXVDSJ5XHgV1V9QFUn4d7KGutnMKq6CxiOmzpmIzBHVdeLyGARGeydtgjYAuR7MQ31M4ZIPfqoK8nq7LPd3cull8J990HLlvDvfyc6KmNMsnj0UXjhhbqB1B3NM5Zzgbki8jfgeGAAbvp8X6nqIlzyCN33WMhnBYaVUUcOkON3bKHmzHE/hyYkrUWmenU3Hczf/uZG7/foAeed58bC1A3mvydjTIqYMwcKC2sHUnc0I++34LqYnsclmdNU9YdAojK+OuEEyM2F225zY16OOgoeftheTTbGBCOS2Y3XeiPc1+BeLz4IaAi85+0zKaBiRfdAf+1a99zl8svdT3u4b4zxWyRdYWcEHoWJmyZN4OWX3W3wVVe5aWEGD4bbb3eDLo0xZm9FMrvxZ6WVeARp/CUCffrAxo1wxRVuxcpmzWD6dFtQzBiz9yLpClvtxznpJifHlVRWrZp7kL96tUssf/87dO4MeXmlTlxgjEkDOTkwaVJuIHVH0hV2VBnPUgSo5lM8JgFat4a33oJZs9xo/SFD2rB6tZs5OR5rzhhj0kskieXICM7JuPeL7rnH/Rw5MrFx+EUEBgxwSyEPHFjAE080YM4cuPVWGDIEykfzYroxJundcw9s3tyAIGaz2etnLF4p8D+05LZwoSvpplo1GDZsMx9+CG3bumcwxxwDr7+e6MiMMX5auBDefbdGIHXHtNCXSX8tWsCrr8K8eW46mK5doVcv2LIl0ZEZY5KdJRZTIhE3oeWGDTBhArzyihtcOWYM/PhjoqMzxiSrqBKLiMQwCbtJdZUru2Ty8cfQty/ccQc0bQrTptnofWPMniJOLCIyDfiviGwVkfdE5HERuTzA2JJalSquZJK6dWHGDHj/fTjiCDfBZdu29vzFmFRUpQpUqhTMX4bR3LGcANRR1QZAL2A+sF8gUaWAf/87c2cLbt8e3n7bTc9fWOiev/Ts6e5ojDGp4d//hjvvXBtI3dEkluXAgQCq+oWqLlLVOwKJyiS9otH7H33kusaWLnUP/K+4Ar79NtHRGWMSqcTEIiJvF9s1FXhDREaKyAkiktGDIm+7zZVMV7kyXHcd5OfDwIHwyCPQuDHcfTfs3Jno6IwxJbntNpg587BA6i4xsajq8cV2PQ3MwQ2qHAosE5HNgUSVApYsccU4tWvD5Mlu9uTjjnMj+I88Ep55xuYfMyYZLVkCq1cfGEjd0XSFFajqzap6h6r2U9UWQPNAojIpq3lzeOkl9x/tQQfBBRe4GZSXLk10ZMaYeIkmseSKyJWhO1T1F5/jMWmiSxe31svMmfD112779NNh3bpER2aMCVo0iaUOMFhEvhSRhSIyXkTOCyowk/r22Qf693dvi911F7zzjpvw8u9/h4KMmwTImMwRzdLE56vqUUAj4CbgY6BDUIEluxo1XDFlq1wZrr0WNm+GESPg6afdgmOjR7vXlY0x8VejBhxwwG+B1B31lC6q+ouqrlbVGap6bRBBpYLnn3fFRK5GDbj3XsjLg3PPdXcxhx/uZlm1N8iMia/nn4dbb10fSN02V5iJu4YN3dovH3wAHTu6u5mmTeHJJ22KGGPSgSWWGI0Z44qJXevWbvTvkiVw8MFwySXQqhW8+CKoJjo6Y9LbmDHw+OONAqnbEkuM3n3XFbP3unSB996DuXPdHcvZZ7slkt94I9GRGZO+3n0X1q8PZpx7NJNQiohcICI3eduHikjGPrw3/hKB3r3d68iPPw5bt0J2NnTrBqtXJzo6Y0w0orljeRToBPTztn8CHvE9IpPRypd3U8Ns2uSmhXn/fTeDcp8+7qG/MSb5RZNYOqrqMGAngKp+D1T0OyAR6SYieSKSLyKjwxwXEXnQO75GRNp4+xuIyFIR2Sgi64sP5jSppUoVGDkSPvkExo51o/lbtHBJZ+vWREdnjClNNInlNxEpByiAiNQCfJ0Fyqv/EaA7brqYfiJSfNqY7kATrwwCJnv7dwHXeGNtjgWGhbnWN/Xru2KCVa2amyxv82YYNsy9Tda4MVx1FXzzTaKjMyZ11a8PtWoFM3lKNInlQdwaLLVFZDzwNjDR53g6APmqukVVfwVmAz2LndMTmKnOcqC6iByiql+p6moAVf0J2AjU8zm+/3n6aVdMfNSpAw884EbxX3ABPPggNGoEN95ogyyNicXTT8MNN2wMpG7RKN7rFJEjga6AAEtU1deoRORcoJuqDvS2++O64IaHnLMQuENV3/a2lwDXqerKkHMaAm8CLVV1j9XZRWQQ7m6HOnXqtJ09e3bYeLZt28b+++/vU+tSRyq0+/PPq/Dkk43IyanN/vv/Rt++W+nVq4AqVfbuJjoV2h6ETG03WNtjbfvJJ5+8SlXbhT2oqhEVYAZQPWT7QGB6pNdH+B3nAdNCtvsDDxU75yXg+JDtJUDbkO39gVVAr0i+s23btlqSpUuXlnjsyitdSUeltTvZfPCB6umnq4Jq7dqq99+vumNH7PWlUtv9lKntVs3ctl95pWrv3ltjvh5YqSX8Xo2mK6yVqhaGJKTvgWOiuD4SBUCDkO36wJeRniMiFYDngWdUdZ7Psf1Jbq4rJrGysmDhQli2DFq2dM9emjSBqVPht2CmQTImLeTmQn5+MHdq0SSWfUTkf6vCiMhBuEW//LQCaCIijUSkItAXWFDsnAXAAO/tsGOBH1T1KxER4Algo6re53NcJsl16vTH4mv168Nll7mFxmbNsmlijIm3aBLLvcC7InKbiNwGLAPu9jMYVd0FDAdexj18n6Oq60VksIgM9k5bBGwB8oHHcatZAhyH6zrrIiK5XunhZ3wm+XXp4u5eXnrJvVE2YIC7k5kzx1ayNCZeIr7jUNWZIrIS6OLt6qWqG/wOSFUX4ZJH6L7HQj4rMCzMdW/jXiowGU4EevRwo/ZfeAFuuskNsGzVCm69Fc46y51jjAlGNFO6VAKygAOAg4Bzi6Z3yURNm7piktc++0CvXvDhh/DMM7Bjh5uHrEMHN/mlTXRpMlnTplC//s+B1B1NV9iLuDEku4DtISUjTZ3qikl+5crBX/8KGzbA9Onw7bfujua44+C11yzBmMw0dSqMHPlxIHVHk1jqq2ofVb1LVe8tKoFEZUwAypeHiy92c45NmeKWRz71VDjpJMjJSXR0xqSPaBLLMhE5OrBIUsygQa6Y1FOxovt3t2kTPPQQ5OfDySdD166wdu0BiQ7PmLgYNAjuuSeY/vxoEsvxwCpvgsg1IrJWRNYEElUK+PhjV0zqqlQJhg9385Ddf7+bsv+KK9rQrZtbH8aYdPbxx1BQsG8gdUeTWIomfzwNOBM4w/tpTEqrUgVGjIAtW2Dw4M2sWgXHHgunnw4rV5Z5uTGmmIgTi6p+Fq4EGZwx8bTfftCnz1Y++QQmToTly6F9e/d68gcfJDo6Y1JHVEsTi8iBItJBRE4sKkEFZkyi7L8/jB7t1oK5/XZ46y1o0wbOOce9umyMKV0041gG4mYMfhm4xfs5Lpiwkl9WlismfR1wANxwA3z6KYwbB0uXun/nvXvD2rUJDs6YvZSVBY0bbwuk7mjuWK4E2gOfqerJuAkoM3appUmTXDHpr1o1uPlml2BuusmNfWnVCs47zz3wNyYVTZoEw4fnB1J3NIllp6ruBDcKX1U/ApoFEpUxSah6dbjlFtdFdsMNsHixSzB9+8L69YmOzpjkEU1iKRCR6sALwGsi8iJ7TmmfMS64wBWTeQ46yD17+fRTGDPGTXh59NEuwWzwffY8Y4JxwQUwfvxRgdQdzVth56hqoaqOA8YC09hz2eCMUVDgislcNWrA+PF/TjAtW0K/frAxmBVfjfFNQQF8802lQOouM7GIyE8i8mNoARYCzwLfBRKVMSmkKMF88glcdx3861/QooWbn+yjjxIdnTHxV2ZiUdWqqnpAsVK1aH88gjQmFdSs6ca/fPopjBoFCxZA8+bwt7+5+cmMyRQ2jsUYn9WsCXfc4e5gRo1ya8I0b+76tC3BmExg41hi1KmTK8aUpFYtl2A+/RRGjoT58y3BmOTRqRO0aPFDIHXbOJYYTZzoijFlqVUL7rzTJZhrrvkjwfTvbwnGJM7EiXDppZ8EUreNYzEmTmrVgrvucl1k11wD8+ZZgjHpKdZxLK9m+jiW3r1dMSZatWuHTzDWRWbiqXdvuOmmFoHUHes4lhuBJ4CzA4kqBXz3nSvGxCo0wVx9tXWRmfj67jv48ccKgdQd1VthRVT1DVVdoKq/+h2QMZmmdm24+27rIjPpI5q3wiqLyNUiMk9EnheRq0SkcpDBGZNJrIvMpIto7lhmAi2Ah4CHgaOAWUEEZUwmsy4yk+qiSSzNVPXvqrrUK4OApkEFluy6dnXFmKBYF5kJUteu0KbN94HUHU1i+UBEji3aEJGOwDt+ByQi3UQkT0TyRWR0mOMiIg96x9eISJtIr/XTjTe6YkzQSuoi698fPv440dGZVHXjjTBgQDCry0cyCeVaEVkDdASWicinIvIJ8C7g65QuIlIOeAToDjQH+olI82KndQeaeGUQMDmKa41JWcW7yObNg6OOggEDLMGY5BLJHcsZwJlAN6ARcBKQ7X0+3ed4OgD5qrrFe+NsNntOzd8TmKnOcqC6iBwS4bW+6d7dFWPiLbSL7OqrYe5cSzAmet27w3XXHR1I3eXLOkFVg7lXCq8esDVkuwB3p1TWOfUivHYPeXl5ZGdnhz1WWFhI9erVwx7LzXU/S7g0pZXW7nSXim0/5hj4/HN45hmYNQvq1IHDDoMqVSKvIxXb7ZdMbXtuLuzatYvs7DLTQNT8r3HvSJh9GuE5kVzrKhAZhOtGo0KFChQWFoYNZvfu3SUe27VrfwAKC7eFPZ7KSmt3ukvVtteoAdWqCV9/XZmvv67If/8rHHjgr9Sps5NKlX4v8/pUbbcfMrXtu3btj6oG0vZkSywFQIOQ7frsOW1MSedUjOBaAFR1KjAVoF27drpy5cqwweTk5JR4N1O0Oycn7OGUVlq70106tP2//3VdZY8+Cps2uXEwY8dC48YlX5MO7Y5VprY9O9vdreXmVo/pepFwf8s70QyQPE9Eqnqfx3oDJduUdV2UVgBNRKSRiFQE+gILip2zABjgvR12LPCDqn4V4bXGpL06deCee2DLFrjySnjuOTjySLj4Yti8OdHRmUwQzevGN6rqTyJyPPAXYAbeG1l+UdVdwHDcWi8bgTmqul5EBovIYO+0RcAWIB94HBha2rV+xhfqjDNcMSZZHXww3Huve8h/xRUwezY0awaXXOKSjslsZ5wBnToFM+FhNF1hu72fpwOTVfVFERnnd0CqugiXPEL3PRbyWYFhkV4blJEj4/Etxuy9gw+G++6Da69168I89hjMnAkXXeS6yBo2THSEJhFGjoScnK3AEb7XHc0dyxciMgU4H1gkIpWivN4Yk0CHHAKTJrm7laFD4emnoUkTGDQI/vOfSokOz6SRaBLD+bhupm6qWggcCFwbRFCpIDs7PV81Numvbl148EH3vOWyy2DGDOjfvyNDhsDWrWVfb9JDdjaMGJEVSN3RJJbTgVdVdZOIjAUeBb4NJCpjTODq1YOHH4b8fOjR4yueeMK9OTZ8OHzxRaKjM6ksqR7eG2Pir0EDuOqqTWzaBBdeCFOmwBFHuDfKvvoq0dGZVBRNYtnj4T1u7IgxJg0cdhhMneqmhbngAnjkETj8cDfx5ddfJzo6k0rs4b0x5k8aNYJp09zU/H37ugf+jRrBddfBt9b5bSKwNw/vDyKDH96ff74rxqSrI46AJ5+EjRvhnHPcaP5Gjdwryt8Hs4yHiaPzz4fs7GBuRSNOLKr6s6rOU9VN3vZXqvpKIFGlgKFDXTEm3TVt6l5NXrcOevSA8ePd2JdbboEffkh0dCZWQ4fC2WeHnfVqr0XVlSUiB4pIBxE5sagEElUK+PlnV4zJFM2bu+lhPvwQunSBcePcHczEibAt/eZiTXs//ww7dwbzNCOaucIGAm/iusNu8X6OCySqFNCjhyvGZJpWrWD+fFi5Ejp3huuvdw/577sPduxIdHQmUj16wOjRrQKpO5p0dSXQHvhMVU8GjgG+CSQqY0zSa9sWFi6Ed9+FrCz39tgRR7i3yX75JdHRmUSKJrHsVNWdACJSSVU/ApoFE5YxJlUceyy88gq88cYfAyybNnVvlv32W6KjM4kQTWIpEJHqwAvAqyLyIiWsd2KMyTwnnuiSyyuvuHnJLr3UPZd55hnYvbvs6036iOatsHNUtVBVxwE3AtMIcE15Y0zqEYFTT3XdYwsWwH77ucGWrVvDvHmgYdd0Nekmmof37URkvoisBh4CJgCrAossyV10kSvGmD2JwJlnwurV7k2yXbugd29o3x4WL7YEkwwuugi6dftPIHVH0xX2DPAk0Bs4EzjD+5mRLLEYU7Z99nED8datc4Mtv/sOuneHk06Ct95KdHSZLVkSyzequkBVP1HVz4pKIFGlgG+/tektjIlU+fLuF1lenntrLD/fPZPp3t3d1Zj4+/Zb+OGHCoHUHU1iuVlEpolIPxHpVVQCiSoFnHuuK8aYyFWs6EZ85+e7KWLef9+9tnzeefDRR4mOLrOcey7cfHOLQOqOJrFcDGQB3XBdYEXdYcYYE5V993VL427ZAjfd5J67tGgBl1wCn3+e6OjM3oomsbRW1XaqeqGqXuyVSwKLzBiT9qpVc3OObdni1n959lm3XPKIETZVfyqLJrEsF5HmgUVijMlYtWq5KWE2bYL+/eGhh9wo/ptvhh9/THR0JlrRJJbjgVwRyRORNSKyVkTWBBWYMSbzNGjgRuyvXw/dusGtt7oEc//9Nk1MKikfxbndAosiBQ0ZkugIjElfRx4J//ynm+hyzBi4+mq34Nitt7oBl+XKJTrC1DdkCKxf/wVQ3fe6oxl5/1m44ntEKaJPH1eMMcFp1w5efdWVWrXcK8tZWW7ySxtkuXf69IEuXYKZR9iWFo7R1q2uGGOCd8opsGKFG8W/c6cb1X/iiW7qGBObrVvh668rBVK3JZYY9e/vijEmPkTcKP4NG2DyZPegv3NnN1VMXl6io0s9/fvDhAlHBVJ30iQWETlIRF4VkU3ezwNLOK+b9wJBvoiMDtl/t4h85L1YMN+bidkYk2YqVIDBg90gy1tvdbMpt2jh9v0nmBlKTJSimYSykoj8VUSuF5GbioqPsYwGlqhqE2CJt108hnLAI0B3oDnQL+QV6FeBlqraCvgYGONjbMaYJLP//nDjjbB5s3sQ/cQTbj2YceNsqeREi+aO5UXcNPm7gO0hxS89gRne5xnA2WHO6QDkq+oWVf0VmO1dh6q+oqq7vPOWA/V9jM0Yk6Rq13bjXjZudMvt3nKLSzCPPeZmVTbxF83rxvVVNchXjuuo6lcAqvqViNQOc049IPSReQHQMcx5lwDPlfRFIjIIGARQp04dcnJywp63bdu2Eo8VFmYBkJOTW9LXpKzS2p3uMrXt6dLuoUMhO7sqjz12BEOGVGfixO0MGrSFzp2/QyT8NenS9mgVFmaxe/fuYNquqhEVYCpwdKTnl1DHa8C6MKUnUFjs3O/DXH8eMC1kuz/wULFzbgDmAxJJTG3bttWSLF26tMRjCxa4ko5Ka3e6y9S2p1u7f/9ddf581aZNVUE1O1t15crw56Zb2yO1YIHq+PFrYr4eWKkl/F6N5o7leOAiEfkE+AUQl5e0VRRJ7JSSjonIf0XkEHV3K4cA4WYKKgAahGzXJ2R5ZBG5EDcxZlev4YE5M2NXojEm+YnA2WfD6afD1KnuuUu7dm5w5YQJboR/pjvzTKha9btA6o7mGUt3oAlwGsEs9LUAuND7fCHumU5xK4AmItJIRCoCfb3rEJFuwHXAWar6s49xhZWXZ684GpPsKlSAYcPcG2SjR7vR/E2buof+P/2U6OgSKy8PPv+8SiB1J9PI+zuAU0VkE3Cqt42I1BWRRV4Mu4DhwMvARmCOqq73rn8YqAq8KiK5IvKYj7Ht4bLLXDHGJL9q1WDiRPfLtFcvuP12N4vyE0/A7t2Jji4xLrsM7ruvWSB1l9kVJiJvq+rxIvITENq9VNQVdoAfgajqd0DXMPu/BHqEbC8CFoU5r7EfcRhj0tdhh8Ezz8AVV8BVV8HAgdC4cVumToWTT050dOmjzDsWVT3e+1lVVQ8IKVX9SirGGBNPHTvCO+/A7Nnw008V6NLF3cls3pzoyNJD0oy8N8aYeBJxEzHOmPE+48e7EfzNm8N119kaMHvLEosxJqNVqvQ7118PH38M/frBXXe5B/xPPgm//57o6FJTRIlFHHtBL8TYsa4YY9JD3brw1FPw/vvQqBFccgkceywsX57oyIIxdiz07x/MyicRJRZvTMgLgUSQok45xRVjTHpp3949f5k1CwoKoFMntw5Muk1wecop0Lbt94HUHe2a9+0DiSIF5ea6YoxJP/vs4wZT5uW58S/PPuu6x+67D377LdHR+SM3F/Lz9w+k7mgSy8nAuyKy2da8hxEjXDHGpK+qVd34l/Xr4YQT4Jpr3AqWS5cmOrK9N2IEPPxwMKM0oh15fwTQhWBG3htjTFJq0sQth7xgAezYAV26uAf9X35Z9rWZKKqR90B1XDI5E6ju88h7Y4xJWiJufq31693cY/PnQ7Nm6dU95pdoFvq6EngGqO2Vp0Xk8qACM8aYZFSlCtx8s0swJ57ousfatoVlyxIdWfKIpivs70BHVb1JVW8CjgUuDSYsY4xJbkcc4brH5s+HwkI47jg3Rcx3wUwYnFKiSSwChE7Xttvbl5EmTHDFGJO5iqbn37ABrr0WZsyAI4+EmTMh2IU79t6ECTBw4JZA6o4msTwJvCci40RkHG753ycCiSoFdO7sijHG7L+/G7G/erVbFvnCC904kU2bEh1ZyTp3hpYtg5m7JuKR98A/gYuB/wO+By5W1UmBRJUCli2zPlVjzJ8dfbQbXDl5Mqxa5bYnTkzOh/vLlsG6dcHMIxzVyHtVXa2qD6rqA6r6QSARpYjrr3fFGGNC7bMPDB4MGzfCGWe43xNt28KKFYmO7M+uvx6mTTs8kLpt5L0xxgTgkENg7lx44QX4v/9z846NHAk/B76+beLZyHtjjAlQz57u1eSBA+Hee6F1a3jjjURHFaxonrEMxkbeG2NM1KpVgylT4PXX3VT82dlw+eWwfXuiIwtGNM9Y7g94zXtjjElrJ58Ma9a4pZEffhhatYI330x0VP6zZywxmjTJFWOMicZ++8EDD7juMBF393LNNbBzZ3zjmDQJhg/PD6TuaJ+xLLdnLE5WlivGGBOLE0+EDz+EIUPcfGNt2rhXlOMlKwsaN94WSN3Rzm58OPaMBYDXXnPFGGNitd9+8Mgj8PLL8OOP7s2xCRNg9+6yr91br70Gq1YdGEjdZSYWERkF/5vduEOx5yuXBRJVCrj9dleMMWZvnXYarF0LvXvDDTe47rFPPw32O2+/HWbNOiyQuiO5Y+kb8nlMsWPdfIzFGGMy1oEHwj/+4ZZEXrPGvZb83HOJjio2kSQWKeFzuG1jjDExEnFLIufmQvPm0LevG/+Saq8lR5JYtITP4baNMcbspUaN3GvI118P06dD+/ZukGWqiCSxtBaRH0XkJ6CV97lo+2i/AhGRg0TkVRHZ5P0M+1RJRLqJSJ6I5IvI6DDHR4qIikhNv2Izxph4q1ABxo+HV191a7x06OC6yVJBmYlFVcup6gGqWlVVy3ufi7Yr+BjLaGCJqjYBlnjbfyIi5YBHcG+oNQf6iUjzkOMNgFOBz32MK6wpU1wxxpggde3qusbat4cBA+DSS/0Z8zJlClx9dd7eVxRGNK8bB60nMMP7PAM4O8w5HYB8Vd2iqr8Cs73ritwPjCIOXXTNmrlijDFBO+QQ93rwmDEwbZobA1NQsHd1NmsGhx66w58AiykfSK2xqaOqXwGo6lciUjvMOfWArSHbBUBHABE5C/hCVT90U5uVTEQGAYMA6tSpQ05OTtjztm3bVuKxZctqANC5c/qtQ1pau9NdprY9U9sNqdX2006DKlVqMnHikbRq9Ts337ye1q1/iKmuZctqsHPnvkCOrzECoKpxK8BrwLowpSdQWOzc78Ncfx4wLWS7P/AQsC/wHlDN2/8pUDOSmNq2baslWbp0aYnHTjrJlXRUWrvTXaa2PVPbrZqabd+wQbVpU9Xy5VWnTo2tjpNOUm3d+vuYYwBWagm/V+N6x6Kqp5R0TET+KyKHqLtbOQT4OsxpBUCDkO36wJe4WZcbAUV3K/WB1SLSQVX/41sDjDEmCRx1FLz/PvTpA4MGQV4e3HknlCuX6MicZHrGsgC40Pt8IfBimHNWAE1EpJGIVMQN3lygqmtVtbaqNlTVhrgE1MaSijEmXVWrBgsXwrBhbp2Xc85JnvEuyZRY7gBOFZFNuDe77gAQkboisghAVXcBw4GXgY3AHFVNobe7jTHGP+XLu+n3H3oIXnoJTjnFvZqcaEnz8F5VvwO6htn/JdAjZHsRsKiMuhr6HZ8xxiSr4cOhbl3461/hhBPcpJYNGpR9XVCSJrGkmlQZqGSMyQy9esHixW4p5M6d3WqVTZqUfP6sWfDuuxuBTr7HkkxdYSmlQYPE/kVgjDHFZWe7BcR27nSrVeaXso5XgwZQu/YvgcRhiSVGzz2XujOPGmPSV1aWu1spSi6bN4c/77nn4PXXawUSgyWWGE2e7IoxxiSbo4+GJUtgxw6XXD77bM9zJk+GBQvqBfL9lliMMSYNtW7tpoH58Uc480z46af4fbclFmOMSVNZWTBnDmzY4N4Yi8eSx2CJxRhj0tppp8GDD7rBlKNGxec77XVjY4xJc0OHwkcfwX33wfHHu1H6QbLEEqO5cxMdgTHGRO6++9yiYTfe6Ma6zJ0L77yzHjjO9++yrrAY1azpijHGpILy5V1SWb8e5s1zv7+qVfstkO+yxBKjp55yxRhjUkWfPm6Br1tvhenTYfHigwP5HkssMbLEYoxJNeXKubuWtWvdjMiWWIwxxuy1vn2hadPwgyb9YonFGGMySNFdy/bt8MMPFQL5DkssxhiTYfr2hRo1YJ99NJD6LbEYY0yGKV8eWraEqlV3BVN/ILVmgEWlLjVmjDHJbdEiePPNNcCJvtdtiSVG++6b6AiMMSZ2++4LlSv/Hkjd1hUWo0cfdcUYY1LRo4/CCy/UDaRuu2OJ0Zw57ufQoYmNwxhjYjFnDhQW1g6kbrtjMcYY4ytLLMYYY3xlicUYY4yvLLEYY4zxlagGM/IyVYjIN0BJs+bUBL6NYzjJIlPbDZnb9kxtN1jbY237YapaK9yBjE8spRGRlaraLtFxxFumthsyt+2Z2m6wtgfRdusKM8YY4ytLLMYYY3xliaV0UxMdQIJkarshc9ueqe0Ga7vv7BmLMcYYX9kdizHGGF9ZYjHGGOMrSyzFiMh5IrJeRH4XkXbFjo0RkXwRyRORvyQqxngQkSwRWS4iuSKyUkQ6JDqmeBKRy71/z+tF5K5ExxNPIjJSRFREaiY6lngRkbtF5CMRWSMi80WkeqJjCpKIdPP++84XkdF+12+JZU/rgF7Am6E7RaQ50BdoAXQDHhWRcvEPL27uAm5R1SzgJm87I4jIyUBPoJWqtgDuSXBIcSMiDYBTgc8THUucvQq0VNVWwMfAmATHExjv99YjQHegOdDP+/3mG0ssxajqRlXNC3OoJzBbVX9R1U+AfCCd/4pX4ADvczXgywTGEm9DgDtU9RcAVf06wfHE0/3AKNy//4yhqq+oatE6vcuB+omMJ2AdgHxV3aKqvwKzcb/ffGOJJXL1gK0h2wXevnQ1ArhbRLbi/mJP27/gwmgKnCAi74nIGyLSPtEBxYOInAV8oaofJjqWBLsE+HeigwhQ4L/LMnKhLxF5DTg4zKEbVPXFki4Lsy+l/6or7Z8D0BW4SlWfF5HzgSeAU+IZX5DKaHt54EDgWKA9MEdEDtc0eDe/jHZfD5wW34jiJ5L/70XkBmAX8Ew8Y4uzwH+XZWRiUdVYfkEWAA1CtuuT4t1Dpf1zEJGZwJXe5j+BaXEJKk7KaPsQYJ6XSN4Xkd9xk/V9E6/4glJSu0XkaKAR8KGIgPvve7WIdFDV/8QxxMCU9f+9iFwInAF0TYc/IkoR+O8y6wqL3AKgr4hUEpFGQBPg/QTHFKQvgZO8z12ATQmMJd5ewLUZEWkKVCTNZ79V1bWqWltVG6pqQ9wvnzbpklTKIiLdgOuAs1T150THE7AVQBMRaSQiFXEvJS3w8wsy8o6lNCJyDvAQUAt4SURyVfUvqrpeROYAG3C3ysNUdXciYw3YpcADIlIe2AkMSnA88TQdmC4i64BfgQvT/C9YAw8DlYBXvTu25ao6OLEhBUNVd4nIcOBloBwwXVXX+/kdNqWLMcYYX1lXmDHGGF9ZYjHGGOMrSyzGGGN8ZYnFGGOMryyxGGOM8ZUlFmOMMb6yxGKMMcZXllhMShORbVGeP05ERgYVT7oTkYYiskNEckP27fbW7VknIv8UkX1Lub6Kd+6vmbTeS6axxGJMlMSJ+v+dWK+LRcDftdlbp6fIDlXNUtWWuJkKShyxrqo7vGtTep49UzpLLCbleX9FbxSRx70VH18RkSohx2/wVst7DWgWsv9q76/sdSIyorT9Id/xKLCaP0/iV3T8IxGZ4a1COFdE9g13nYhcICLve3+5TxGRciKyn4i8JCIfet/bJ9y+kO9aF/LdI707sYi+K8w/v9ki8py3TMBnInL6XvzreAto7NX7gois8v6dZNKUQBnPEotJF02AR7wVHwuB3gAi0hY3yd4xuJVB24fsvxjoiJse/1IROaak/d53NANmquoxqvpZmBiaAVO9VQh/BIYWvw7YF+gDHOf95b4b+BtuVdIvVbW195f/4hL2lSWS7yquNbBFVTt6x2+O4Hv24M0r1x1Y6+26RFXbAu2AK0SkRiz1mtRjicWki09UNdf7vApo6H0+AZivqj+r6o/8MYvr8d7+7aq6DZjnnVvSfoDPVHV5KTFsVdV3vM9Pe3UVv64r0BZY4T2n6AocjvtlfIqI3CkiJ6jqDyXsK0sk3/U/3p1dTeAWb9cG3Fo0Rcdvi+A7q3j1r8QtafyEt/8KEfkQtyJjA1zyNxnAZjc26eKXkM+7gSoh2+FmWg232FFp+wG2lxFD8e8p2g69ToAZqrrHipze3VIPYKKIvKKqt4bbh5tdO/SPwsolxFjid4VoCWxS1Z3edhvgQy+eg4nsd0TRc5PQtmTjFobrpKo/i0hOsThNGrM7FpPu3gTO8d5GqgqcGbL/bO85yH7AObjnAyXtj8ShItLJ+9wPeDvMOUuAc0WkNoCIHCQih4lIXeBnVX0atxR0m3D7vDr+C9QWkRoiUgm3OFU4Yb+r2Dmtvbgre+29BbfuPbjuw9wI215cNeB7L6kcietWNBnC7lhMWlPV1SLyHO4X5Gd4ScLb/xR/LNY2TVU/AAi3X0QaRvB1G4ELRWQKbmG0yUDtYvFsEJGxwCveW1u/AcNwv4jvFrda5W/AEODoMPtQ1d9E5FbgPeAT4KMS2l7Sd4U+H2qNW4Y3BzgAmBDSnZeF6wqMxWJgsIisAfJw3WEmQ9h6LMb4wEs8C72H7ClDRN4ELlXVvDDHnvCO/R6yryE+tFNEPgXaqWpar8yZqawrzJjMdgQlLDutqn8PTSqe3UA1CRkgGQ2vSzIXqAAUr9ukCbtjMcYY4yu7YzHGGOMrSyzGGGN8ZYnFGGOMryyxGGOM8ZUlFmOMMb6yxGKMMcZXlliMMcb46v8B2j8V85WwKqcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot error in mass balance\n",
    "pi = np.linspace(pi_min, pi_max, 100)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.plot(pi, f(pi), 'b')\n",
    "ax.axvline(x=pi_min, color='b', **{'linestyle': 'dashed'})\n",
    "ax.axvline(x=pi_max, color='b', **{'linestyle': 'dashed'})\n",
    "ax.axhline(color='k')\n",
    "\n",
    "ax.set_ylabel(r'Error in mass balance $\\dot{m}$ [kg/s]')\n",
    "ax.set_xlabel(r'Indoor pressure $p_i$ [Pa]')\n",
    "ax.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c20ac20-e752-4e87-bff2-4ad722797f68",
   "metadata": {},
   "source": [
    "### Numerical solution\n",
    "\n",
    "Alternativelly, the solution to the equation\n",
    "\n",
    "$$f(p_i) = 0$$\n",
    "\n",
    "can be found numerically."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "033906e8-4963-433a-ac81-2fd62c8af01c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "pi = -4.45 Pa\n"
     ]
    }
   ],
   "source": [
    "# numerical solution\n",
    "root = float(fsolve(f, np.mean([pi_min, pi_max])))\n",
    "print(f\"pi = {root:.2f} Pa\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f88d82d2-244f-42e6-af9a-e68ff6ea540d",
   "metadata": {},
   "source": [
    "### Check the mass balance\n",
    "\n",
    "By knowing the indoor pressure $p_i$, the mass flows  $\\dot m_{ni}$, $\\dot m_{si}$ and $\\dot m_{ic}$  can be found and the mass balance\n",
    "\n",
    "$$ \\dot m_{ni} + \\dot m_{si} - \\dot m_{ic} = 0$$\n",
    "\n",
    "checked."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "99073457-1203-4d7b-864b-cff69fc912eb",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Verify mass balance:\n",
      "mass flow: north: 0.019 kg/s, south: 0.019 kg/s\n",
      "mass flow chimney: 0.038 kg/s\n",
      "error in mass flow balance: 0.000 kg/s\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# verification of mass balance, kg/s\n",
    "mn = ρo * Ko * (Pdn - root)**no\n",
    "ms = ρo * Ko * (Pds - root)**no\n",
    "mc = ρi * Kc * (root + Psi - Pdc - Pso)**nc\n",
    "print(\"Verify mass balance:\")\n",
    "print(f\"mass flow: north: {2 * mn:.3f} kg/s, south: {2 * ms:.3f} kg/s\")\n",
    "print(f\"mass flow chimney: {mc:.3f} kg/s\")\n",
    "print(f\"error in mass flow balance: {2 * mn + 2 * ms - mc:.3f} kg/s\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48938e42-383d-4a2b-a3bd-93ebb5052198",
   "metadata": {},
   "source": [
    "## Considering the wind effect\n",
    "\n",
    "### Dynamic pressure\n",
    "\n",
    "The dynamic pressures are:\n",
    "- $P_{dn} = C_pn \\frac{1}{2} \\rho_o v^2$ - north face\n",
    "- $P_{ds} = C_ps \\frac{1}{2} \\rho_o v^2$ - south face\n",
    "- $P_{dc} = C_pc \\frac{1}{2} \\rho_o v^2$ - chimney"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "a83dd661-e4a0-4f4e-b5fe-dfc0abce9d43",
   "metadata": {},
   "outputs": [],
   "source": [
    "v = 18 * 1000 / 3600            # m/s, wind speed\n",
    "\n",
    "# dynamic pressures, Pa\n",
    "Pdn = Cpn * 0.5 * ρo * v**2     # north face\n",
    "Pds = Cps * 0.5 * ρo * v**2     # south face\n",
    "Pdc = Cpc * 0.5 * ρo * v**2     # chimney"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4c422902-9fec-4f5f-90ba-1a87281589d5",
   "metadata": {},
   "source": [
    "### Domain of variation of indoor pressure\n",
    "The domain of variation of $p_i$ in order to have the flows shown in Figure 1 is:\n",
    "\n",
    "$$P_{si} - (P_{dc} + P{sc}) \\le p_i \\le \\min (P_{dn},P_{ds})$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "aac14797-e1c6-43f1-9957-c1aba346426b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Pressure domain with wind: -21.23 Pa < pi < -9.04 Pa\n"
     ]
    }
   ],
   "source": [
    "pi_min, pi_max = -(Psi - Pdc - Pso), min(Pdn, Pds)\n",
    "print(f\"Pressure domain with wind: {pi_min:.2f} Pa < pi < {pi_max:.2f} Pa\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9a3d88c0-b5c1-4326-99c5-7f32db3256f1",
   "metadata": {},
   "source": [
    "### Graphical solution of nolinear equation for indoor pressure\n",
    "\n",
    "The solution of this non-linear equation can be found graphically."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "60584731-ffe1-43e9-adce-eec79721c660",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZYAAAEKCAYAAAAxXHOuAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAA5EklEQVR4nO3deZyN5f/H8dfH2JciW77Inn3LWMsu+74Touxrv1YllVJKWpA1JFEMZd+KDMkSMlkjlAj11TdFyvr5/XGO73dilnNmzj33OTOf5+NxP5x7O/f7Mua63Nt1iapijDHGBEoqtwMYY4xJXqxhMcYYE1DWsBhjjAkoa1iMMcYElDUsxhhjAiq12wHcliNHDi1YsKDf+x369RAAxbMXD3Ai3/z5559kypTJlWMHkpUjuFg5gouT5UhsHbZr166zqpozpnUpvmEpWLAgO3fu9Hu/OrPrABDZMzKwgXwUGRlJnTp1XDl2IFk5gouVI7g4WY7E1mEicjy2dXYpzBhjTEAF3RmLiDQGxgNhwAxVffWm9SWA94B7gBGqOs7XfQPp2VrPOvXVxhjjOCfrsKBqWEQkDJgE3A+cBHaIyDJVPRBts/8AQ4HWCdg3YBoUbuDE1xpjTJJwsg4LtkthVYAjqnpMVS8D84FW0TdQ1V9UdQdwxd99AynqTBRRZ6Kc+npjjHGUk3VYUJ2xAHmBE9HmTwJVA72viPQF+gLkzp2byMhIv4M+EvUIAG9XeNvvfQPhwoULCcodbKwcwcXKEVycLIeTdViwNSwSwzJfe8n0eV9VnQ5MBwgPD9eEPHWR9YesAK49eWJPvQQXK0dwsXLEz8k6LNguhZ0E8kebzwecSoJ9jTHGBEiwNSw7gGIiUkhE0gKdgWVJsK/fzl86z9HfjmLDDhhjzD8F1aUwVb0qIoOBtXgeGZ6lqvtFpL93/VQRuRPYCdwGXBeRR4BSqvpHTPs6lfXClQuc/OMk+37ZR9ncZZ06jDHGhJygalgAVHUVsOqmZVOjfT6D5zKXT/s6ZWyDsbRe0Jo1R9ZYw2KMCTmv1H/Fse8OtkthIaNViVaUzVWW1UdWux3FGGP8ViN/DWrkr+HId1vDkkBbTmyhbK6ybP5xM+cvnXc7jjHG+GXLiS1sObHFke+2hiWBnln/DPv/vZ8r16+w4YcNbscxxhi/PLP+GZ5Z/4wj320NSyLcnu52MqXJxJoja9yOYowxQcMalkQQEeoXrs/qI6vtsWNjjPGyhiWRmhRtwg/nfuDwr4fdjmKMMUHBGpZEalSkEYBdDjPGGK+ge48lVLzd+G0ACmUrRPHsxVl9ZDXDqg1zN5QxxvjoRh3mBGtYEqjCnRX++7lJ0SZM2TmFPy//Saa0oT/OtjEm+YtehwWaXQpLoHXH1rHu2DoAOpbuyKVrlxizeYzLqYwxxjfR67BAszOWBBq9aTTgGYWtev7q9Cjfg7FfjqVLmS6UzlXa5XTGGBO36HVYoNkZS4CMu38cWdJlof/K/lzX627HMcYY11jDEiA5M+Vk3P3j2PzjZmbtnuV2HGOMcY01LAHUs0JPahWoxROfPcGp8zbGmDEmZbKGJYBEhGnNp3H52mU6LerElWtX3I5kjDFJzm7eJ9C05tNiXF4iRwlmtJhB10+68tS6p3iz0ZtJnMwYY+IXWx0WCNawJFDxHMVjXdelbBe2ntzKW9veolq+anQs3TEJkxljTPziqsMSyy6FJdDyQ8tZfmh5rOvHNRxH9XzVeWjpQ+z9eW8SJjPGmPjFV4clhjUsCfTG1jd4Y+sbsa5PG5aWhR0Wcnv622kyrwknfj+RhOmMMSZu8dVhiWENi4Py3paX1Q+s5vzl8zSZ14Tf/vrN7UjGGOM4a1gcVi53OZZ0WsLhXw/Tan4r/r76t9uRjDHGUdawJIG6heoyp80cvvjxC9pFtOPS1UtuRzLGGMfE+1SYiNzhw/dcV9VziY+TfHUu05k/Lv1BvxX96LCwA4s6LiJtWFq3YxljTMD58rjxKe8kcWwTBtwVkEQh4oM2H/i9T99Kfbl2/RoDVw2k06JORLSPIE1YGgfSGWNM3BJSh/nKl4bloKpWjGsDEdkdoDwhI//t+RO034DKA7h6/SpD1wylXUQ7IjpEkD51+gCnM8aYuCW0DvOFL/dYqgdom2Rlwb4FLNi3IEH7Dqk6hElNJ7H88HKafdiMC5cvBDidMcbELTF1WHzibVhU9W8AEekgIlm8n0eKyCcick/0bQJBRBqLyCEROSIiw2NYLyIywbt+z40M3nX/JyL7RWSfiHwkIo6dCkzZOYUpO6ckeP+BlQcyp/UcNv6wkQZzGvCfv/4TwHTGGBO3xNZhcfHnqbCRqnpeRO4DGgLvAwFNJSJhwCSgCVAK6CIipW7arAlQzDv1vZFBRPICQ4FwVS2D575P50DmC7Tu5buzsMNCdp/ZzX2z7uPH3390O5IxxiSaPw3LNe+fzYApqroUCPRjTVWAI6p6TFUvA/OBVjdt0wqYox7bgKwikse7LjWQQURSAxnxPHQQ1NqUbMPabms5df4U1WZU45sz37gdyRhjEsWfTih/EpFpQAPgNRFJR+Dfg8kLRO/75CRQ1Ydt8qrqThEZB/wI/AV8qqqfxnQQEemL52yH3LlzExkZ6XfQc+fOASRo35i8WeZNhu8bTo0ZNRhVahThd4THuf2FCxcCdmw3WTmCi5UjuDhZjkDXYf+gqnFOeG7MC54zgLZAMe/yPEDD+Pb3ZwI6ADOizXcHJt60zUrgvmjz64FKQDbgcyAnkAZYAnSL75iVKlXShKj9Xm2t/V7tBO0bmxO/n9Cyk8tq2KgwnbJjSpzbbtiwIaDHdouVI7hYOYKLk+VIbB0G7NRY6lVfzlgexHPf4zCwBjjvbZBOA6cT1ard6iQQ/Rm4fNx6OSu2bRoA36vqvwFE5BOgBjA3wBkBWNRxUcC/M99t+fjyoS/p8nEXBqwcwLdnv+WNhm8Qlios4McyxqRsTtRhN/jyVFh/Vb0HeAHPWcFsEdkqIq+ISC3vDfdA2QEUE5FCIpIWz833ZTdtswzo4X06rBrwu7eR+xGoJiIZRUSA+sDBAGb7hxwZc5AjY46Af2+WdFlY2nkp/1ft/xi/fTzNPmxmnVcaYwLOqToM/LhHoqrfqupbqtoYqAdsxnPpanugwqjqVWAwsBZPoxChqvtFpL+I9Pdutgo4BhwB3gUGevfdDiwCvgb24inb9EBlu9nsqNnMjprtyHeHpQrjzUZvMr35dD7//nMqv1uZ/b/sd+RYxpiUyck6LEEjSKrqX3gq+FWBjQOqesv3qurUaJ8VGBTLvs8Dzwc6U0xu/EB6Vujp2DH6VOpD6VylaRfRjmozq/F+6/dpW7KtY8czxqQcTtZhPjcsIvJoDIt/B3apalTAEpl/qJG/Bjv77KRdRDvaRbTjyRpP8nL9l92OZYwxsfLnceFwoD+ex33z4nlctw7wrog8Gfho5oa8t+VlY8+NDAgfwNgtYz1v6l+2N/WNMcHJn4YlO3CPqj6mqo/haWhyArWAng5kM9GkS52Oyc0mM6f1HL766Sv67OrDhu83uB3LGGNu4U/DchdwOdr8FaCA936LjVyVRLqX78723tvJFJaJBh804KWNL3Ht+rX4dzTGmCTiz837D4FtIrLUO98C+EhEMgEHAp4syK16IODPLfisbO6yTKs0jXm/z+O5yOfYeHwjH7T5gDxZ8sS/szHG4Gwd5s8ZyyqgD3AOz037/sBuVf1TVR9wIFtQy5gmIxnTZHTt+BnCMvBBmw+Y0WIGW05sofzU8qz+brVreYwxocXJOsyfhuVd4LKqjlfVt/H0LvysI6lCwOQdk5m8Y7KrGUSEh+95mJ19d3Jn5jtp+mFTHl37KJeu2pVJY0zcnKzD/GlY2gPvi0hJEemD58XEho6kCgER+yOI2B/hdgwASuUsxVd9vmJQ5UG8te0tqsyoYi9UGmPi5GQd5s+b98fwdLHyMZ5GpqGq/u5IKuO39KnT807Td1jeZTmnz58m/N1wJm6fyHW97nY0Y0wKE2/DIiJ7vSM17sHTZcodQEFgu3eZCSLN727O3gF7qVuwLkPXDKXx3Mb89MdPbscyxqQgvjwV1tzxFCagcmfOzcquK5m2axqPffoYZaaUYXLTyXQp28XtaMaYFMCX3o2PxzUlRUjjPxGhf3h/ovpFUSJHCbp+0pWOCzty9uJZt6MZY5I5Xy6FfR2IbZKbyJ6RRPaMdDtGvIplL8YXvb7glXqvsOTbJZSeXJol3y5xO5YxxmVO1mG+3LwveeMeSyzTXsCZTv1NQKROlZqnaz7Nrr67yJslL20WtKHbJ9349eKvbkczxiRDvtxjKeHDNimuT5FxW8YB8HiNx11O4ruyucuyvfd2XvniFUZ/MZp1x9YxpdkU2pRs43Y0Y0wSc7IOS/Q9Fu90MuDJgtyKwytYcXiF2zH8liYsDc/XeZ6dfXaSJ0se2ka0pdOiTvx84We3oxljkpCTdZg/L0iaZKT8neX5qvdXvFT3JZZ8u4RSk0sxd89cPOOoGWNMwlnDkoKlCUvDs7WeZXe/3dyd/W66L+5O0w+bcvycPexnjEk4vxoWEUnjVBDjnlI5S7G512bGNx7PF8e/oPTk0ozfNt664zfGJIjPDYuIzAB+FpETIrJdRN4VkSEOZgtqGdJkIEOaDG7HCJiwVGEMrTqU/QP3U6tALR5Z+wjVZ1Yn6kyU29GMMQ5wsg7zZzyWmkBuVb0iInmB8kA5R1KFgNUPJM8u6gtkLcDKriuZv28+j6x9hPDp4fxftf/jhTovkCltJrfjGWMCxMk6zJ9LYduAbACq+pOqrlLVV52JZdwkInQp24WDgw7Ss0JPxm0dR+nJpVl5eKXb0YwxISDWhkVENt+0aDqwUUQeF5GaInK7s9GC20sbX+KljS+5HcNRd2S4gxktZ7Cx50YypslI84+a0z6ivXVqaUwy4GQdFmvDoqr33bRoLhCB5/LZQGCLiBx1JFUIWP/9etZ/v97tGEmiVoFaRPWP4uV6L7Pyu5WUmFSCt7a+xdXrV92OZoxJICfrMH8uhZ1U1edV9VVV7aKqpYFSjqQyQSdtWFqeqfkM+wfup+ZdNXn000cJnx7O1hNb3Y5mjAky/jQsUSIyLPoCVbUxcFOYwtkKs7LrShZ1WMTZi2epMasGDy992HpNNsb8lz8NS26gv4icEpEVIvKyiHQIdCARaSwih0TkiIgMj2G9iMgE7/o9InJPtHVZRWSRiHwrIgdFpHqg8xnPzf12pdpxcNBBHq/+OHP2zOHuiXczbec0e/fFGOPX0MQdVbUkUAh4DjgMVAlkGBEJAyYBTfBcZusiIjdfbmsCFPNOfYEp0daNB9aoagk8j0MfDGS+6LJnzE72jNmd+vqQkCVdFl5v+DpR/aIof2d5+q/sT9UZVdl+crvb0Ywx8XCyDvPnPRbgv5e/vvZOgVYFOKKqxwBEZD7QCjgQbZtWwBz1dGq1zXuWkgf4E6gF9PTmvAxcdiAjAB93/Niprw45pXOV5vMenzN/33we+/Qxqs2sxsMVH2ZM/THkzJTT7XjGmBg4WYdJMHU6KCLtgcaq2ts73x2oqqqDo22zAnhVVTd759cDTwFX8TwSfQDP2couYJiq/hnDcfriOdshd+7clebPn+9ouZxw4cIFMmfO7HaMW1y8epE5x+ew6KdFZAjLQK+CvWj1r1aESViM2wdrOfxl5QguVg7n1a1bd5eqhse4UlWDZgI6ADOizXcHJt60zUrgvmjz64FKQDiexqWqd/l44KX4jlmpUiVNiOGfDdfhnw1P0L6BsGHDBteO7YsDvxzQ++fcr7yAlp1cVj8/9nmM2wV7OXxl5QguVo74JbYOA3ZqLPVqsPVufBLIH20+H3DKx21O4nkk+sYF/kXAPThk68mtbD1pj9rGpmTOkqzttpZPOn7C+cvnqTenHh0XdrSek40JEk7WYf50Qiki0k1EnvPO3yUiAb15D+wAiolIIRFJC3QGlt20zTKghzdPNeB3VT2tqmeAEyJS3Ltdff55b8YkMRGhTck2HBh4gBfrvMiKwysoMakEz294notXLrodzxjjEH/OWCYD1YEu3vnzeJ7gChhVvQoMBtbieaIrQlX3i0h/Eenv3WwVcAw4AryLpxeAG4YA80RkD1ABeCWQ+UzCZEiTgZG1R3Jo8CFal2jNi5tepPg7xZm/b74NLGZMMuTPU2FVVfUeEdkNoKq/ec8qAkpVV+FpPKIvmxrtswKDYtk3Cs+9FhOE8t+en4/afcTA8IEMWzOMLh93ocxtZXiv+HuE/8t+bMYkF/6csVzxvmeiACKSE7juSKoQkO+2fOS7LZ/bMUJSzQI12dFnB++2eJef/vqJyu9WptfSXpw6f/PtNGOMU5ysw/w5Y5kALAZyicjLQHtgpCOpQsDctnPdjhDSwlKF0fue3uT5Tx426kbGbx/Pwv0LGX7fcB6r/liyGkTNmGDkZB3mz5v384AngTHAaaC1qkY4FcykDJlSZ2Ls/WM5MPAAjYo2YuSGkRR/pzgf7v3Q7r8YE6L8eSrsfeCMqk5S1XeAMyIyy7lowe2RNY/wyJpH3I6RbBS5owgfd/yYDQ9uIEfGHDzwyQNUn1mdLSe2uB3NmGTJyTrMn3ss5VT13I0ZVf0NqBjwRCEi6kyUjQfvgDoF67Cjzw7ea/UeP/7+I/fOupdOizrx/W/fux3NmGTFyTrMn4YllYhkuzEjIneQgL7GjIlPWKowelboyeEhh3mu1nMsP7ScEpNK8ORnT3Lu73NuxzPGxMOfhuUNYKuIvCQiLwFbgNediWUMZE6bmVF1R/HdkO/oWrYr47aMo+iEokzcPpEr1664Hc8YEwt/bt7PAdoCP3untt5lxjgq7215ea/Ve+zqu4vyd5Zn6JqhlJlShsUHF9sNfmOCkD8379PheZv9NuAOoP2N7l1Soruz383d2e92O0aKUjFPRdZ1X8eKLisIkzDaRrSl1uxaNv6LMQngZB3mzz2SpcDveLqjT/FDEk9vMd3tCCmSiNDs7mY0KtqIWbtn8dyG56g2sxodSnVgTP0xFLmjiNsRjQkJTtZh/jQs+VS1sWNJjPFD6lSp6Vup73/vvby+5XWWfLuEAeEDGFl7JDky5nA7ojEplj8377eISFnHkoSYvsv70nd5X7djpHiZ02bmhTovcGTIEXpV6MU7O96hyIQijPlijPWgbEwcnKzD/GlY7gN2icghEdkjInu9vQinSId/PczhXw+7HcN45cmSh2ktprF3wF7qFKzDM58/w90T72bm1zO5dv2a2/GMCTpO1mH+NCxNgGJAQ6AF0Nz7pzFBo1TOUiztvJRNPTeR77Z89F7em3JTy7H80HJ7gsyYJOLP48bHY5qcDGdMQtUsUJOtD29lUYdFXLl2hZbzW1J7dm22nrBRP41xml9DE4tINhGpIiK1bkxOBTMmsUSEdqXasX/gfiY3nczhXw9TY1YN2i5oy7dnv3U7njHJlj/vsfQGNuEZ3XGU988XnIkV/CrcWYEKd1ZwO4bxQZqwNAyoPIAjQ4/wUt2XWHdsHWUml6Hv8r789MdPbsczxhVO1mH+PG48DKgMbFPVuiJSAk8DkyK93fhttyMYP2VOm5lnaz1Lv0r9ePmLl5m8YzIf7PmAYVWH8dS9T5EtQ7b4v8SYZMLJOsyfS2F/q+rf4HkLX1W/BYo7E8sY5+TMlJO3G7/N4SGHaV+qPWO/HEvhCYV5bfNr9oiyMQHgT8NyUkSyAkuAdSKyFEixY8l2+6Qb3T7p5nYMkwgFsxbkgzYfENU/ihr5azB8/XCKTSzG9F3TrZNLk+w5WYf581RYG1U9p6ovAM8CM4BWjqQKASf/OMnJP066HcMEQLnc5VjZdSUbe26kYNaC9FvRj9KTSxOxP4Lret3teMY4wsk6LN6GRUTOi8gf0SdgBfAh8KsjqYxxQa0CtdjcazPLOi8jXep0dFrUifDp4aw9stbegTHGD/E2LKqaRVVvu2nKcmN5UoQ0JqmICC2KtyCqXxQftPmA3/7+jcbzGlP3/br2DowxPrL3WIyJQViqMLqV68ahwYeY2GQiB88epMasGrT8qCV7f97rdjxjgpq9x5JA1fNVp3q+6m7HMA5LG5aWwVUGc3ToUUbXHc2m45soP7U83T7pxtH/HHU7njEJ5mQdZu+xJNCYBmPcjmCSUOa0mRlRawQDKg9g7JdjmbB9Agv2L6B3xd48W+tZt+MZ4zcn67Cge49FRBp7e1A+IiLDY1gvIjLBu36PiNxz0/owEdktIisCnc2YOzLcwasNXuXo0KP0q9SPmbtnUnRiUaYencrZi2fdjmdMUEjoeyyfOfEei4iEAZPw9KRcCugiIqVu2uxGL8vFgL7AlJvWDwMOBjJXTNpFtKNdRDunD2OCVJ4seXin6TscGnyITqU7sfDkQgqPL8wLkS/wx6U/3I5nTLycrMMS+h7LSGAm0DrAeaoAR1T1mKpeBuZz67syrYA56rENyCoieQBEJB/QDM87No769eKv/HrRnrZO6QplK8Ts1rOZGT6T+4vcz6iNoyg0vhCvf/m6vcVvgpqTdZg/91j+S1U3BjqIV17gRLT5k0BVH7bJC5wG3gaeBLLEdRAR6YvnbIfcuXMTGRnpd9Bz584BJGjfQLhw4YJrxw6k5FKOHJqDIbmG0DBDQ2Z+P5Mn1z3Ja5teo9td3WiWpxlpUqVxO6JPksvPw8oRPyfrMJ8bFhFJDwzEM5KkApuBKTfuuwSIxLDs5jfTYtxGRJoDv6jqLhGpE9dBVHU6MB0gPDxc69SJc/MYZf0hKwAJ2TcQIiMjXTt2ICW3ctShDv3ox6bjmxjx+QjGHxnP0rNLeb7283Qr143UqRL0f7kkk9x+HqHOyXI4WYf5c49lDlAamAi8A5QEPghwnpNA/mjz+bj1Pk5s29wLtBSRH/BcQqsnInMDnM8Yn9QqUItNPTex5oE1ZM+QnV5Le1F2SlnrJsakCP40LMVV9WFV3eCd+gJ3BzjPDqCYiBQSkbRAZ2DZTdssA3p4nw6rBvyuqqdV9WlVzaeqBb37fa6qjvUSWb9QfeoXqu/U15tkQERoVLQRO/rs4JOOnxAmYXRa1Il7pt1jQyUb1zlZh/lzXr5bRKp5b5gjIlWBLwMZRlWvishgPC9fhgGzVHW/iPT3rp8KrAKaAkeAi0CvQGbw1cjaI904rAlBIkKbkm1oWbwl8/fN5/nI52k5vyVV81ZldL3R1C9UH5GYrvAa4xwn67B4GxYR2YvnPkcaPGcKP3rnCwAHAh1IVVfhaTyiL5sa7bMCg+L5jkggMtDZjEmMsFRhPFDuATqW7sicb+YwauMo7v/gfmoXqM3oeqO576773I5oTED4cimsOdACaAwUAmoDdbyfmzmWLMg1mdeEJvOauB3DhKA0YWl4+J6H+W7Id0xsMpFDvx6i5ns1aTy3MTt+2uF2PJNCOFmH+dK78fG4JkdShYC/rvzFX1f+cjuGCWHpUqf7bz9krzV4jZ2ndlJlRhVaz2/Nnp/3uB3PJHNO1mF+9W5sjAm8jGky8uS9T3Js2DFerPMiG37YQPmp5em0qBMH/+14JxLGBJw1LMYEidvS3cbI2iP5YdgPjKg5gpWHV1JmShl6LO5hPSmbkOJPt/kdRCSL9/OzIvLJzR1AGmMSL1uGbIyuN5rvh33P/1X7PxYeWEjxd4rTZ1kffvz9R7fjGRMvf85YRqrqeRG5D2gEvM+tHUCmGM3vbk7zu5u7HcMkYzkz5WRcw3EcG3qMAeEDmLNnDsUmFmPwqsGcOh/Q/l9NCuRkHeZPw3LN+2czPF25LAXSBj5SaHi8xuM8XuNxt2OYFCBPljxMbDqR74Z8R8/yPZm2axpFJhThsbWP8cufv7gdz4QoJ+swfxqWn0RkGtARWCUi6fzc3xiTCHfdfhfTWkz7b1f9b29/m8LjC/P0uqetp20TVPxpGDrieSO+saqeA7IBTzgRKhTUmV2HOrPruB3DpECFsxVmduvZHBh4gJbFW/Lal69RaHwhnt/wPOf+Pud2PBMinKzD/GlYmgGfqep3IvIsMBmwIfOMcUnxHMX5sN2H7Bmwh4ZFGvLiphcpNL4QozeN5vyl827HMymY3bw3JsSVyVWGRR0X8XXfr6l5V01GbhhJofGFGPvlWP68/Kfb8UwKZDfvjUkmKuapyLIuy9jeezuV81bmqXVPUXhCYd7a+pb1EmGSlN28NyaZqZK3CqsfWM3mXpspk6sMj376KEUnFmXSV5O4dPWS2/FMCpCYm/d3kIJv3ncs3ZGOpTu6HcOYWN17172s77GeDQ9uoHC2wgxePZhiE4sxfdd0rly74nY84zIn6zCfx2NR1YvAJ9HmT+MZZz5FGlh5oNsRjPFJnYJ12NRzE+uOrWPkhpH0W9GPVze/ynO1nwuJ4ZKNM5ysw/y6lCUi2USkiojUujE5FSzYXbxykYtXLrodwxifiAj3F7mfrQ9vZWXXlWTLkI1eS3tRalIp5u2Zx7Xr1+L/EpOsOFmH+dNXWG9gE57LYaO8f77gSKoQ0HReU5rOa+p2DGP8IiI0LdaUnX12srjTYtKnTk+3xd0oN7UcC/cv5LpedzuiSSJO1mH+nLEMAyoDx1W1LlAR+LcjqYwxjhIRWpdoTVT/KBa0X4Cq0nFRR/rs6sPig4vxDNRqTML407D8rap/A4hIOlX9FijuTCxjTFJIJanoWLojewfsZV7beVy+fpm2EW0JfzeclYdXWgNjEsSfhuWkiGQFlgCfichSwLpYNSYZCEsVRteyXZldeTbvtXqP3/76jeYfNaf6zOp8dvQza2CMX3xuWFS1jaqeU9UXgJHADKCVU8GMMUkvTMLoWaEnhwYfYlrzaZw6f4qGcxtSe3ZtIn+IdDueCRE+P2coIuHACKCAdz8BXgbKORMtuPWs0NPtCMY4Jk1YGvpW6suD5R9kxtczePmLl6n7fl3qFarHS3Vfokb+Gm5HNInkZB3mzwPs8/C8ELkXSPGPjljDYlKCdKnTMajKIB6q+BDTdk1jzOYx3DvrXhoVacSLdV+kSt4qbkc0CeRkHebPPZZ/q+oyVf1eVY/fmBxLFuTOXjzL2YvWubNJGTKkycAj1R7h2NBjvNbgNXae2knVGVVp+VFLos5EuR3PJICTdZg/DcvzIjJDRLqISNsbkyOpQkD7iPa0j2jvdgxjklSmtJl48t4n+X7Y97xU9yW++PELKk6rSPuI9uz/Zb/b8YwfnKzD/GlYegEVgMZAC+9kg74bkwJlSZeFZ2s9y/fDvue5Ws/x6dFPKTulLF0/7srhXw+7Hc+4zJ+Gpbyqhqvqg6rayzs9FOhAItJYRA6JyBERGR7DehGRCd71e0TkHu/y/CKyQUQOish+ERkW6GzGmH/Kmj4ro+qO4vth3/PUvU+x9NBSSk4qSc8lPTn22zG34xmX+NOwbBORUo4lAUQkDJgENAFKAV1iOGYToJh36sv/Bhu7CjymqiWBasAgp/MaYzyyZ8zOmAZjODb0GMOqDmPB/gUUf6c4/Zb348TvJ9yOZ5KYPw3LfUCU92xij4jsFZE9Ac5TBTiiqsdU9TIwn1vflWkFzFGPbUBWEcmjqqdV9WsAVT0PHATyBjifMSYOuTPn5s1Gb3J06FH6VerHe1HvUXRiUYauHsrp8ym2M/QUR3x9o1ZECsS0PJBPholIezzjvfT2zncHqqrq4GjbrABeVdXN3vn1wFOqujPaNgXxdJhZRlX/iOE4ffGc7ZA7d+5K8+fP9zvr5798DkC9XPX83jcQLly4QObMmV05diBZOYJLoMtx5u8zzD0+l9VnVpMmVRpa/6s1nfN3JmvarAE7Rkzs5xG/xNZhdevW3aWq4TGuVNWgmYAOwIxo892BiTdtsxK4L9r8eqBStPnMwC6grS/HrFSpkoaiDRs2uB0hIKwcwcWpchz59Yj2WNxDU41KpZlfyawj1o/Q/1z8jyPHUrWfR1IAdmos9WqwDS18EsgfbT4ft/ZHFus2IpIG+BiYp6qf4KATv5+wa8fG+KjIHUV4v/X77Buwj2bFmvHyFy9TaHwhRm8azflL592OlyI5WYcFW8OyAygmIoVEJC3QGVh20zbLgB7ep8OqAb+r6mkREWAmcFBV33Q6aPfF3em+uLvThzEmWSmZsyTz288nql8UdQrWYeSGkRQaX4hxW8bx15W/3I6XojhZhwVVw6KqV4HBeAYROwhEqOp+EekvIv29m60CjgFHgHeBG+Nr3ovn0lk9EYnyTjYSlzFBqPyd5VnSeQnbe2+n0r8q8cRnT1BkQhEmfTWJS1cvuR3PJJI/nVCmA9oBBaPvp6ovBjKQqq7C03hEXzY12mcFBsWw32Y8HWMaY0JElbxVWNttLZuOb+LZz59l8OrBvL7ldZ6r/Rw9yvcgdSp/ujM0wcKfM5aleB71vQr8GW0yxphEqVWgFht7bmRtt7XkypSLh5c9TKlJpfho70c2XHII8qdhyaeqnVR1rKq+cWNyLJkxJkURERoWacj23ttZ0mkJ6VOnp+snXSk/tTxLvl1ig42FEH/OM7eISFlV3etYmhDyWPXH3I5gTLIkIrQq0YoWxVsQsT+C5zY8R5sFbaiStwqj646mQeEGeJ7VMYnhZB3m75v3uxx+8z5ktCjeghbFW7gdw5hkK5WkonOZzhwYdIAZLWZw+vxpGs5tSL059dhyYovb8UKek3WYPw3LjT66GvK/no1TbM166OwhDp095HYMY5K91KlS8/A9D/PdkO8Y33g8B/59gHtn3UvzD5vbWDCJ4GQd5s+Y98djmhxJFQL6rehHvxX93I5hTIqRLnU6hlYdyrGhxxhTfwxbTmyh4rSKdPm4i3XVnwBO1mHxNiwicqNPrvMi8ke06byI3NIPlzHGOClT2kwMv284x4YdY0TNESw/tJxSk0rRZ1kf6w0jSMTbsKjqfd4/s6jqbdGmLKp6m/MRjTHmVlnTZ2V0vdEcHXqUQZUHMWfPHIpNLMajax/l3OVzbsdL0YLqzXtjjPFX7sy5Gd9kPIcHH6Zr2a6M3z6erl915YXIF/jjkl1UcYM1LMaYZKFA1gLMajWLfQP2UTlbZUZtHEWRCUV4a+tb/H31b7fjpSg+vcfi7eAxn6raBUyvZ2s963YEY0wMSuYsyajSo8hULBPPfP4Mj376KG9te4tRdUbRvXx36ybGy8k6zKczFm//XEscSxGCGhRuQIPCDdyOYYyJReW8lfms+2es676OOzPfyUPLHqLclHIsPrjY3uLH2TrM3zHvKzuSIgRFnYmyZ+iNCQH1C9dne+/tfNzxY67rddpGtKXGrBps/GGj29Fc5WQd5k/DUhfYKiJH7c17eGTNIzyy5hG3YxhjfCAitC3Zln0D9zGjxQxO/H6COu/XodmHzdjzc8qsxpysw/x9874IUA97894YE4Kiv8U/tsFYtpzYQoWpFeixuAc/nPvB7XjJhl9v3gNZ8TQmLYCsKfnNe2NM6MqQJgNP3PsEx4Ye44kaT7DwwEKKv1OcR9c+ytmLZ92OF/J8blhEZBgwD8jlneaKyBCnghljjNOyZcjGa/e/xuHBh+lWthvjt4+nyIQijPliDBevXHQ7Xsjy51LYw0BVVX1OVZ8DqgF9nIlljDFJJ//t+ZnZaiZ7+u+hdoHaPPP5M9w98W7e2/0e165fczteyPHngW4Bov8NXyMFDwX8Sv1X3I5gjAmw0rlKs6zLMjYd38QTnz3BQ8se4s1tbzK2wVgaF22crMaBcbIO8+eM5T1gu4i8ICIvANuAmY6kCgE18tegRv4abscwxjigVoFabHt4GxHtI/jryl80/bApDT5owNenv3Y7WsA4WYf51LB437xfCPQC/gP8BvRS1bcdSRUCtpzYYoMNGZOMiQgdSnfgwKADTGg8gW/OfEOl6ZXosbgHP/7+o9vxEs3JOsyvN+9V9WtVnaCq41V1tyOJQsQz65/hmfXPuB3DGOOwtGFpGVJ1CEeHHuWpe58iYn8Ed0+8m6fXPR3SnVw6WYfZm/fGGOOD29PfzqsNXuXQ4EO0L9WeV798laITijJlxxSuXr/qdrygYm/eG2OMHwpkLcDctnPZ0WcHJXOWZOCqgZSbUo7V3612O1rQ8OceS3/szXtjjAEg/F/hRD4YyeJOi7ly/QpNP2xKo7mN2PfLPrejuc6feyxv2Zj3xhjzPyJC6xKt2T9wP281eouvfvqK8lPLM3DlQP7957/djueaoLvHIiKNReSQiBwRkeExrBcRmeBdv0dE7vF130B6u/HbvN34bScPYYwJEWnD0vJItUc4MuQIgyoPYvqu6RSbWIw3t77J5WuX3Y4XIyfrMH/vsWxz8h6LiIQBk/B0eFkK6CIipW7arAlQzDv1Bab4sW/AVLizAhXurODU1xtjQlD2jNmZ0GQCewfspXr+6jz26WOUnVKWVd+tcjvaLZysw/x5876JIwn+qQpwRFWPAYjIfKAVcCDaNq2AOd7Lc9tEJKuI5AEK+rDvLQ4dOkSdOnX8Dvrb378BkC19Nr/3DYRz586RNWtWV44dSFaO4GLlCKyyf5XlyH+O0GxcM7JlyEbRO4qSMU1Gn/d3shxO1mHxNiwi8qSqjlXV4yLSQVUXRlv3ChDIB6HzAtGHPz4JVPVhm7w+7guAiPTFc7ZDmjRpOHfunN9Bj144CkDRzEX93jcQrl27lqDcwcbKEVysHIGVilQUy1SMs5fOcubvM+w4tYOcaXOSO31uwiQs3v2dLIeTdZgvZyydgbHez0/jeQP/hsYEtmGJqSOem8cQjW0bX/b1LFSdDkwHCA8P1507d/qTEYA6s+sAENkz0u99AyEyMjJBZ1rBxsoRXKwczvn5ws+M+HwEs3bPIlWmVIy9fyzdynUjlcR+R8LJciS2Dour3zRf7rFILJ9jmk+sk0D+aPP5gFM+buPLvsYY44rcmXMzo+UMtvfeTsGsBXlwyYPUfK8mu08nv05MfGlYNJbPMc0n1g6gmIgUEpG0eM6Wlt20zTKgh/fpsGrA76p62sd9jTHGVZXzVmbLw1uY1XIW3/36HeHvhjNk1RDO/X3O7WgB40vDUl5E/hCR80A57+cb82UDGUZVrwKDgbXAQSBCVfeLSH8R6e/dbBVwDDgCvAsMjGvfQOYzxphASCWp6FWxF4eHHGZg+EAm75xM8XeKM+ebOXieSwpt8d5jUdX47zAFkKquwtN4RF82NdpnBQb5uq9TpjWflhSHMcYkY1nTZ2Vi04k8VPEhBq4ayINLHmTW7llMaTaFkjlLOnpsJ+swf95jMdEUz1Gc4jmKux3DGJMMVMxTkS8f+pLpzaez5+c9lJ9anhHrR3Dp2iXHjulkHWYNSwItP7Sc5YeWux3DGJNMpJJU9KnUh28Hf0uXsl14ZfMrPLTzIdYfW+/I8Zysw6xhSaA3tr7BG1vfcDuGMSaZyZUpF++3fp/1PdYjIjT4oAEPLnmQXy/+GtDjOFmHWcNijDFBqF6heswKn8WImiP4cO+HlJxUkoj9ESFxc98aFmOMCVJpU6VldL3R7Oq7iwJZC9BpUSfaLGjD6fOn3Y4WJ2tYjDEmyJXLXY6tD2/l9ftfZ+3RtZSeXJp5e+YF7dmLNSzGGBMCUqdKzeM1HieqXxQlcpSg2+JutFnQhp8v/Ox2tFv407uxieaDNh+4HcEYkwIVz1GcL3p9wdvb3mbE5yMoO6UsM1rOoGXxln59j5N1mJ2xJFD+2/OT//b88W9ojDEBFpYqjMdqPMauvrvIe1teWs1vRe9lvblw+YLP3+FkHWYNSwIt2LeABfsWuB3DGJOClc5Vmu29tzP83uHM2j2L8OnhfHPmG5/2dbIOs4YlgabsnMKUnVPcjmGMSeHShqVlTIMxrOuxjj8u/UHVGVWZvGNyvDf2nazDrGExxphkoF6henzT/xvqFarHoFWD6La4GxevXHQlizUsxhiTTOTMlJMVXVfwcr2X+WjvR1SfWZ1jvx1L8hzWsBhjTDKSSlLxTM1nWPXAKk78foJK0ysR+UNk0mZI0qMZY4xJEo2LNmZn353kyZyHRnMbsXD/wvh3ChB7jyWBFnVc5HYEY4yJU+Fshdn80GZaftSSTos6cebCGYZUHQI4W4dZw5JAOTLmcDuCMcbE644Md/BZ98/o+klXhq4ZSrrU6ehbqa+jdZg1LAk0O2o2AD0r9HQ1hzHGxCdDmgws6rCI5h81Z9iaYdyb/152nNoBOFOHWcOSQNawGGNCSViqMGa3mk35qeXp/HFnsqXPRipJ5UgdZjfvjTEmhcidOTfvt36ffb/s4+h/jjp2HGtYjDEmBWlUtBGPV3+cUxdOcfbiWUeOYQ2LMcakMC/Xf5nsGbKTOpUzd0PsHosxxqQwacPSUiZXGce+3xqWBFr1wCq3IxhjTII5WYdZw5JAGdNkdDuCMcYkmJN1mN1jSaDJOyYzecdkt2MYY0yCOFmHBU3DIiJ3iMhnIvKd989ssWzXWEQOicgRERkebfnrIvKtiOwRkcUiktXJvBH7I4jYH+HkIYwxxjFO1mFB07AAw4H1qloMWO+d/wcRCQMmAU2AUkAXESnlXf0ZUEZVywGHgaeTJLUxxph/CKaGpRXwvvfz+0DrGLapAhxR1WOqehmY790PVf1UVa96t9sG5HM2rjHGmJgE08373Kp6GkBVT4tIrhi2yQuciDZ/Eqgaw3YPAbEO5iwifYG+ALlz5yYyMtLvsOfOnQNI0L6BcOHCBdeOHUhWjuBi5QguTpbDyTosSRsWEVkH3BnDqhG+fkUMy/4xsLOIjACuAvNi+xJVnQ5MBwgPD9c6der4ePj/yfpDVgASsm8gREZGunbsQLJyBBcrR3BxshxO1mFJ2rCoaoPY1onIzyKSx3u2kgf4JYbNTgL5o83nA05F+44HgeZAfVVVfLBr166zInLcpwLElLtXTG1dksgBONMfQ9KycgQXK0dwcbwciajDCsS2IpguhS0DHgRe9f65NIZtdgDFRKQQ8BPQGegKnqfFgKeA2qp60deDqmrOROZ2hYjsVNVwt3MklpUjuFg5gkuoliOYbt6/CtwvIt8B93vnEZF/icgqAO/N+cHAWuAgEKGq+737vwNkAT4TkSgRmZrUBTDGGBNEZyyq+itQP4blp4Cm0eZXAbf0RaCqRR0NaIwxxifBdMZi/DPd7QABYuUILlaO4BKS5RAf73EbY4wxPrEzFmOMMQFlDYsxxpiAsoYlxMTV2aaIPO3tnPOQiDRyMWacRKSDiOwXkesiEh5teRoReV9E9orIQREJ6v7eYiuHd105EdnqXb9XRNK7lTM+cZXDu/4uEbkgIo+7kc9Xcfy7ul9Ednl/DrtEpJ6bOeMTz7+rkPgdt4Yl9MTY2aa3M87OQGmgMTDZ22lnMNoHtAU23bS8A5BOVcsClYB+IlIwibP5I8ZyiEhqYC7QX1VLA3WAK0meznex/TxueAtYnXRxEiy2cpwFWnj/XT0IfJDUwfwU27+rkPkdD5rHjY1vVPXTaLPbgPbez62A+ap6CfheRI7g6bRzaxJHjJeqHgQQueWNXwUyeSvmDMBl4I+kTee7OMrRENijqt94t/s1iaP5JY5yICKtgWPAn0mbyn+xlUNVd0eb3Q+kF5F03t+VoBPHzyNkfsftjCW0PcT//icZUwedeZM8UeIswlOBnQZ+BMap6n/cjZQgdwMqImtF5GsRedLtQAkhIpnw9GYxyu0sAdQO2B2sjUo8QuZ33M5YglBcnXWq6lLvNjd3thlvB51JyZcyxKAKcA34F5AN+EJE1qnqMYdixiuB5UgN3AdUBi4C60Vkl6qudyhmvBJYjlHAW6p6IaazGTcksBw39i0NvIbnjNJVCSxHUP2Ox8UaliAUV2edEGtnm3F20JnU4itDLLoCa1T1CvCLiHwJhOO5FOOKBJbjJLBRVc8CeLskugfPAHauSGA5qgLtRWQskBW4LiJ/q+o7AQ3nhwSWAxHJBywGeqjq0cCm8l8i/l0Fze94XOxSWIiJ1tlmy5s621wGdBaRdN5OOosBX7mRMRF+BOqJRyagGvCty5kSYi1QTkQyeu8X1QYOuJzJb6paU1ULqmpB4G3gFTcblYTyPjm5EnhaVb90OU5ihMzvuDUsoSfGzja9nXFG4KnA1gCDVPWaezFjJyJtROQkUB1YKSJrvasmAZnxPBWzA3hPVfe4FDNesZVDVX8D3sRThijga1Vd6VrQeMTx8wgpcZRjMFAUGOn9nYmSmAcSDApx/LsKnd9x69LFGGNMINkZizHGmICyhsUYY0xAWcNijDEmoKxhMcYYE1DWsBhjjAkoa1iMMcYElDUsxhhjAsoaFhPSROSCn9u/EOzjigQzESkoIn+JSFS0Zde8Lx3uE5GFIpIxjv0zeLe9LCI5kiS0SXLWsBjjJ2+XM37/7iR0v4Rw+FhHVbVCtPm/VLWCqpbBM9RB/9h2VNW/vPsGZR9XJjCsYTEhz/u/6IMi8q535L1PRSRDtPUjvCPurQOKR1v+qPd/2ftE5JG4lkc7xmTga/7ZGeCN9d+KZwTMPSKyyNtX2C37iUg3EfnK+z/3aSISJiKZRGSliHzjPW6nmJZFO9a+aMd+3Hsm5tOxYvj7my8iC0Rku4gcF5FmifhxfIGn+xREZIl4RmzcLyJ9E/GdJsRYw2KSi2LAJO+IjefwjLuBiFTCM+peRTyj8lWOtrwXnh58qwF9RKRibMu9xygOzFHViqp6PIYMxYHp3tE9/wAG3rwfkBHoBNzr/Z/7NeABPCMCnlLV8t7/+a+JZVl8fDnWzcoDx1S1qnf98z4c5xbi6XCzCbDXu+ghVa2Ep4fqoSKSPSHfa0KPNSwmufheVaO8n3cBBb2fawKLVfWiqv6Bp4dY8IyXslhV/1TVC8An3m1jWw5wXFW3xZHhRLTec+d6v+vm/erjGXZ5h/c+RX2gMJ7KuIGIvCYiNVX191iWxceXY/2X98wuB/8bzOsAnrFwbqx/yYdjZvB+/048PVTP9C4fKiLf4BnpND+ext+kADYei0kuoo8IeA3P0MY3xNTTamwjV8U1olV8w/PefJwb89H3E+B9VX36lgN7zpaaAmNE5FNVfTGmZXgGeIv+n8L0sWSM9VjRlAG+U9W/vfP3AN9489yJb3XEjfsm0ctSB2gAVFfViyISeVNOk4zZGYtJ7jYBbbxPI2UBWkRb3tp7HyQT0AbP/YHYlvviLhGp7v3cBdgcwzbr8QyelQtARO4QkQIi8i/goqrOBcYB98S0zPsdPwO5RCS7iKTDM+hbTGI81k3blPfmTu8t7yjgLe+6ini6/U+I24HfvI1KCTyXFU0KYWcsJllT1a9FZAGeCvI43kbCu3w2/xsoaYaq7gaIabmIFPThcAeBB0VkGvAdMAX4x7gfqnpARJ4FPvU+tXUFGISnIn5dRK57lw0AysawDFW9IiIvAtuB74llMLQ4jhX9/lB5PMNbRwK34RnM68blvAp4LgUmxBqgv4jsAQ7huRxmUggbj8WYAPA2PCu8N9lDhohsAvqo6qEY1s30rrsebVlBAlBOEfkBCL8xfLNJXuxSmDEpWxE8Z1e3UNWHozcqXteA2yXaC5L+8F6SjALSADd/t0km7IzFGGNMQNkZizHGmICyhsUYY0xAWcNijDEmoKxhMcYYE1DWsBhjjAkoa1iMMcYElDUsxhhjAur/ATXzy/cOxRJXAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot error in mass balance\n",
    "pi = np.linspace(pi_min, pi_max, 100)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.plot(pi, f(pi), 'g')\n",
    "ax.axvline(x=pi_min, color='g', **{'linestyle': 'dashed'})\n",
    "ax.axvline(x=pi_max, color='g', **{'linestyle': 'dashed'})\n",
    "ax.axhline(color='k')\n",
    "\n",
    "ax.set_ylabel(r'Error in mass balance $\\dot{m}$ [kg/s]')\n",
    "ax.set_xlabel(r'Indoor pressure $p_i$ [Pa]')\n",
    "ax.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f2812f53-e590-48fc-973c-b865e3ab902a",
   "metadata": {},
   "source": [
    "### Numerical solution\n",
    "\n",
    "Alternativelly, the solution to the equation\n",
    "\n",
    "$$f(p_i) = 0$$\n",
    "\n",
    "can be found numerically."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "76fa348a-e59e-4f66-989b-daa72fad8de3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "pi = -10.50 Pa\n"
     ]
    }
   ],
   "source": [
    "# numerical solution\n",
    "root = float(fsolve(f, np.mean([pi_min, pi_max])))\n",
    "print(f\"pi = {root:.2f} Pa\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9f286a8a-e7d5-40ab-b599-0c2d27ba8e3c",
   "metadata": {},
   "source": [
    "### Check the mass balance\n",
    "\n",
    "By knowing the indoor pressure $p_i$, the mass flows  $\\dot m_{ni}$, $\\dot m_{si}$ and $\\dot m_{ic}$  can be found and the mass balance\n",
    "$$ \\dot m_{ni} + \\dot m_{si} - \\dot m_{ic} = 0$$\n",
    "checked."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "9c553c84-997c-4fb9-ad00-0937474ec8ee",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Verify mass balance:\n",
      "mass flow: north: 0.045 kg/s, south: 0.009 kg/s\n",
      "mass flow chimney: 0.054 kg/s\n",
      "error in mass flow balance: 0.000 kg/s\n",
      "\n"
     ]
    }
   ],
   "source": [
    "mn = ρo * Ko * (Pdn - root)**no\n",
    "ms = ρo * Ko * (Pds - root)**no\n",
    "mc = ρi * Kc * (root + Psi - Pdc - Pso)**nc\n",
    "print(\"Verify mass balance:\")\n",
    "print(f\"mass flow: north: {2 * mn:.3f} kg/s, south: {2 * ms:.3f} kg/s\")\n",
    "print(f\"mass flow chimney: {mc:.3f} kg/s\")\n",
    "print(f\"error in mass flow balance: {2 * mn + 2 * ms - mc:.3f} kg/s\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "41a0b71a-a98d-40c1-b0bd-7336a4d8e1b3",
   "metadata": {},
   "source": [
    "## References\n",
    "B. Rock, S.J. Emmerich, S. Taylor (2017) *F16 Ventilation and Infiltration* in [ASHRAE Fundamentals Handbook](https://ia802804.us.archive.org/23/items/ashraehandbook2017hvacfundamentalssi/2017%2C%20ashrae%20fundamental%20handbook%20SI.pdf) ISBN 10: 1939200598  ISBN 13: 978-1939200594"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0a0e8140-5cd1-4e28-a9a1-23023ac77e86",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "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.7"
  },
  "toc-autonumbering": true
 },
 "nbformat": 4,
 "nbformat_minor": 5
}