{ "metadata": { "name": "" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "raw", "metadata": {}, "source": [ "Content provided under a Creative Commons Attribution license, CC-BY 4.0; code under MIT license. (c)2014 Lorena A. Barba, Olivier Mesnard. Thanks: NSF for support via CAREER award #1149784." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "[@LorenaABarba](https://twitter.com/LorenaABarba)" ] }, { "cell_type": "heading", "level": 5, "metadata": {}, "source": [ "Version 0.3 -- April 2014" ] }, { "cell_type": "heading", "level": 1, "metadata": {}, "source": [ "Source-vortex panel method" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In [Lesson 9](http://nbviewer.ipython.org/urls/github.com/barbagroup/AeroPython/blob/master/lessons/09_Lesson09_flowOverCylinder.ipynb) of _AeroPython_, you learned to use a _source panel method_ to represent a circular cylinder, and in [Lesson 10](http://nbviewer.ipython.org/urls/github.com/barbagroup/AeroPython/blob/master/lessons/10_Lesson10_sourcePanelMethod.ipynb) we used it for a symmetric airfoil at zero angle of attack. But what if we want the airfoil to generate some lift? If we place the airfoil at a non-zero angle of attack, we _should_ get lift, but will a source-panel representation be able to give you lift? Remember the _Kutta-Joukowski theorem_?\n", "\n", "\n", "Historically, the first panel method ever developed was a source-sheet method. At the time, Douglas Aircraft Company was concerned with calculating the flow around bodies of revolution, and it was only later that the method was extended to lifting surfaces. (See the reference below for a nice historical account.)\n", "\n", "A *source-panel method* leads to a solution with no circulation, therefore no lift. The objective of this lesson is to start with the source panel method we implemented in the previous lesson and add some *circulation* so that we may have a lift force. We introduce an important concept: the **Kutta-condition** that allows us to determine what the right amount of circulation should be." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Reference\n", "\n", "* Smith, A.M.O., The Panel Method: Its Original Development. In _Applied Computational Aerodynamics_, Vol. 125, edited by P.A. Henne, published by AIAA (1990). [Read it on Google Books.](http://books.google.com/books?id=5Ov2tHj0wxoC&lpg=PA3&ots=SnUiqcdEnb&dq=The%20Panel%20Method%3A%20Its%20Original%20Development&pg=PA3#v=onepage&q&f=false)" ] }, { "cell_type": "heading", "level": 2, "metadata": {}, "source": [ "A lifting-body panel method" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If we were to simply increase the angle of attack in the freestream and calculate the flow with a source sheet only, the rear stagnation point will not be located at the trailing edge. Instead, the flow will bend around the trailing edge and the stagnation point will be somewhere on the top surface of the airfoil. This is not a physically possible solution.\n", "\n", "For example, using the source-sheet panel method of [Lesson 10](http://nbviewer.ipython.org/github/barbagroup/AeroPython/blob/master/lessons/10_Lesson10_sourcePanelMethod.ipynb) with an angle of attack $\\alpha=4^\\circ$ (using 40 panels), and plotting the streamlines in an area close to the trailing edge, we get the following plot:\n", "\n", "
![image](files/resources/StreamlinesTrailingEdge.png)
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As you can see, the streamlines behave strangely at the trailing edge. We know experimentally that the flow leaves the trailing edge of an airfoil smoothly, so this must be wrong. Well, it's just that we can't use purely sources to calculate potential flow of an airfoil at non-zero angle of attack\u2014we need circulation. But how do we obtain circulation?\n", "\n", "### The Kutta condition\n", "\n", "The *Kutta-condition* states that the pressure below and above the airfoil trailing edge must be equal so that the flow does not bend around it, but leaves tangentially. The rear stagnation point must be exactly at the trailing edge.\n", "\n", "It's natural to be a little perplexed by this. How can we justify this seemingly arbitrary condition? Remember that potential-flow theory completely ignores fluid viscosity, so if we are leaving out this physical effect, we shouldn't be surprised that the theory needs some adjustment for those situations when viscosity does play a role. A real viscous fluid is not able to turn around a sharp corner like an airfoil trailing edge without separating there. The Kutta condition allows us to correct potential-flow theory so that it gives a solution closer to reality.\n", "\n", "Remember [Lesson 6](http://nbviewer.ipython.org/github/barbagroup/AeroPython/blob/master/lessons/06_Lesson06_vortexLift.ipynb), where we studied lift on a cylinder by combining a doublet and a freestream, plus a vortex. That's when we learned that **lift always requires circulation**. If you experimented with the circulation of the point vortex (which you *did*, right?), you found that the stagnation points moved along the cylinder.\n", "\n", "Like for the circular cylinder, the amount of circulation we add to an airfoil will move the stagnation points along the surface. And if we add just the right amount, the rear stagnation point can be made to coincide with the trailing edge. This amount of circulation makes the flow a physically relevant solution. And this amount gives the correct lift!\n", "\n", "To implement the Kutta-condition in our panel method we need to add one more equation to the system, giving the circulation that moves the stagnation point to the trailing edge. By placing a vortex-sheet with the same constant strength at every panel, we can add the circulation to the flow with just one more unknown.\n", "\n", "How do we enforce this in our code? We can re-use most of the code from Lesson 10, and enforce the Kutta-condition while adding circulation to the flow. Previously, we discretized the geometry into `N` panels, with a constant source strength on each one (varying from panel to panel), and applied a Neumann boundary condition of flow tangency at the `N` panel centers. This led to a linear system of `N` equations and `N` unknowns that we solved with the NumPy function `linalg.solve()`. In the lifting-body case, we will instead have `N+1` equations and `N+1` unknowns. Read on to find out how!\n" ] }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ "Discretization into panels" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's get the preliminaries out of the way. We need to import our favorite libraries, and the function `integrate` from SciPy, as in Lesson 10." ] }, { "cell_type": "code", "collapsed": false, "input": [ "from math import *\n", "import numpy as np\n", "from scipy import integrate\n", "import matplotlib.pyplot as plt" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 1 }, { "cell_type": "markdown", "metadata": {}, "source": [ "We start by importing the NACA0012 geometry from a data file, and we plot the airfoil:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "# reads of the geometry from a data file\n", "with open ('./resources/naca0012.dat') as file_name:\n", " x, y = np.loadtxt(file_name, dtype=float, delimiter='\\t', unpack=True)\n", "\n", "# plots the geometry\n", "%matplotlib inline\n", "\n", "val_x, val_y = 0.1, 0.2\n", "x_min, x_max = x.min(), x.max()\n", "y_min, y_max = y.min(), y.max()\n", "x_start, x_end = x_min-val_x*(x_max-x_min), x_max+val_x*(x_max-x_min)\n", "y_start, y_end = y_min-val_y*(y_max-y_min), y_max+val_y*(y_max-y_min)\n", "\n", "size = 10\n", "plt.figure(figsize=(size, (y_end-y_start)/(x_end-x_start)*size))\n", "plt.grid(True)\n", "plt.xlabel('x', fontsize=16)\n", "plt.ylabel('y', fontsize=16)\n", "plt.xlim(x_start, x_end)\n", "plt.ylim(y_start, y_end)\n", "plt.plot(x, y, color='k', linestyle='-', linewidth=2);" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAnEAAAB+CAYAAABLREfEAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xt0E2X+P/B3mvSW9N6kadP0QinScgqlKwocF211UfCC\n62W1q6i4yAJnXYVlF1zdr4rHG15WXXvOgscLK0cRd3UBBaoLpV3AlopQLwVckDa9p9e0JaGXpPP7\ng1/mJPQWCkwy6ft1znMyz2Rm+syH0Hz6PM/MKARBEEBEREREshLg7QYQERER0fljEkdEREQkQ0zi\niIiIiGSISRwRERGRDDGJIyIiIpIhJnFEREREMqTydgOklpubi5KSEm83g4iIiGhU11xzDYqLi4d8\nTzHe7hOnUCggx1NetGgRNm7c6O1mjCuMufQYc+kx5tJjzKUn55iPlLdwOJWIiIhIhpjEyURqaqq3\nmzDuMObSY8ylx5hLjzGXnr/GnEmcTOTm5nq7CeMOYy49xlx6jLn0GHPp+WvMmcQRERERyRCTOCIi\nIiIZ4tWpRERERD6KV6cSERER+Zlxd7NfuSouLvbbiZm+SuqYC4IAh8MBu90+bLnYvcgKhQIqlQqB\ngYHiq7MolUooFIqL+vNGw8+59Bhz6THm0vPXmDOJIxrBwMAATp8+ja6uLnR3d7u9Wq1W2Gw28dVZ\nXOtnzpxBT08P+vr60NvbO2zp7++Hw+Hw9ukO4prYnZvkDbXOWQ8ODkZoaKjHJSQkBKGhoTh27Bii\noqKG3EapVHo7HEREPsXn5sQVFhZixYoVcDgceOihh7BmzZpB2zzyyCPYtWsX1Go1Nm7ciJycHABn\n7wMTEREBpVKJwMBAlJeXD9qXc+LGH0EQcPr0aTQ3N6OlpQVtbW1ob29HR0cHOjo6xGXX166uLnR1\ndeH06dOStlWpVEKlUg1ZlEolAgIu7gyIgYEBsZevv7/frQwMDFzUn3WhAgMDodFoEB4ejrCwMISF\nhYnL5756sk6j0UCl4t+xROTbRspbfOo3mMPhwMMPP4zdu3cjMTERV1xxBRYsWIDMzExxm507d+Lk\nyZM4ceIEDh48iOXLl6OsrAzA2RMtLi5GTEyMt06BJGS1WtHQ0OBWGhsb0dzcLCZsLS0taG5uRm9v\n75h/jvNLPyIiwu01LCwMarUaGo0GarVaLK710NBQBAcHj1qCgoK8Mnw5EmeC50zqhkr0zl3nrPf0\n9Ii9kGfOnBm1eLJdf38/LBYLLBbLRTvHkJAQhIeHIzw8HFFRUYiMjBRfXZeHe42MjERgYOBFaw8R\n0fnwqSSuvLwc6enp4p2V8/PzsW3bNrckbvv27XjggQcAADNnzoTFYoHZbIZerwcAv+1l89fx/OH0\n9vaitrYWJpMJJpMJNTU14qszYevq6vL4eGq1GnFxcdDpdIiNjUVMTAxiYmIQHR2N6Ohocdn5GhkZ\niW+//Rbz5s0bt8N4AQEBCAoKQlBQkGQ/c7jPuSAI6Ovrg9VqxenTp9Hd3Y3Tp0+7LY+0brj3enp6\n0NPTg5aWljG3Wa1WD5vgRUVFiZ+1oYparfZ64j7efrf4AsZcev4ac59K4urr65GUlCTWjUYjDh48\nOOo29fX10Ov1UCgU+MUvfgGlUomlS5diyZIlkrWdzl97eztOnjw5qFRVVaGpqWnU/YODg2EwGNxK\nQkIC9Ho9dDqdmLTpdDpoNJrzbt/JkyfHbQLnaxQKhdhrebF62gVBwJkzZ8Q5jp2dnejs7ITFYhnx\n9dx1zvmPjY2N592GoKAgMaFz/eNitBIeHu715I+IvM+nkjhPfykN19u2f/9+GAwGtLS0YO7cucjI\nyMCcOXMGbbdo0SKxty8qKgrTp08XM/Ti4mIAYP0i1YuKitDU1ISIiAgcPXoURUVFMJlMMJvN6Ojo\nwHCUSiWMRiMiIiKg1+sxc+ZMJCcnw2KxQKfT4dZbb0V0dDRKSkpGbU91dfWY2p+bm+v1+I23unOd\nFD9PoVCI82Zd34+NjcUdd9zh0fH27t2Lnp4eTJ06FZ2dnSgqKoLVakVycjI6Oztx+PBhdHd3Iyws\nDO3t7Th16hS6urrQ19eHtrY29PT0oKmpyaM/WlwplUrExsYiNDQUUVFRSE9Ph1arhdVqRWRkJGbP\nng2dTofq6mpERkZiwYIFCA0NHfZ8nLz9788665eqniuj3+fO5erqaozGpy5sKCsrw9NPP43CwkIA\nwAsvvICAgAC3ixuWLVuG3Nxc5OfnAwAyMjJQUlIiDqc6rV27FmFhYVi1apXbel7YcOl0dHSgoqIC\nR44cQUVFBSorK3Hs2DGcOXNmyO3DwsKQnp4+qEyYMAEGg4GTzsnvnTlzBu3t7eddxnLBjVqthlar\nhU6ng1arFYtr3dlzHRcXh+jo6It+IQ0RnT/ZXNgwY8YMnDhxAtXV1TAYDNiyZQs2b97sts2CBQtQ\nUFCA/Px8lJWVISoqCnq9HjabDQ6HA+Hh4bBarfjyyy/x1FNPeelMLr5iHxvP7+joQHl5OQ4ePIgj\nR47gyJEjMJlMQ25rMBgwZcoUsWRmZmLy5MmIi4vz6SEhX4v5eDDeYh4aGorExEQkJiae137OnrzW\n1la0traipaVlxOWWlhbYbDbU1NSgpqbGo5+hUqnEhC4uLg56vX7IZWcJCQkZSwjGpfH2OfcF/hpz\nn0riVCoVCgoKcMMNN8DhcGDx4sXIzMzEhg0bAABLly7FjTfeiJ07dyI9PR0ajQbvvfceAKCpqQm3\n3347AMBut+Pee+/F9ddf77Vz8ScDAwOorKxEaWkpysrKUFpaiuPHjw/aLiQkBNOmTUNOTg5ycnIw\ndepUTJkyBVFRUV5oNZH/CgoKQkJCAhISEjza3nmbnaESvW+++QZqtVqsO6/utlgsaGxs9HiuX0RE\nxIjJXnx8POLj45GQkICwsLALOX0i+v98ajhVChxOHZ0gCDh69Cj27t2LvXv3oqSkBG1tbW7bBAcH\n4/LLL8fMmTNx+eWXIycnB5dddhmHQIn8RG9vr3iLnubmZpjN5hGX7Xa7x8fWaDRISEgQkzrXV9dl\nnU7Hi4to3Bspb2ESRwDOXin6xRdf4PPPP8eePXtgNpvd3k9KSsJVV12F2bNnY9asWcjOzkZwcLCX\nWktEvkQQBPF2T8Mles4LOBobG9HT0+PRcQMCAhAXF+dRwjeWK9CJ5IBJnAu5JnGXYjz/xx9/xNat\nW7Fjxw4cOHDA7Q798fHxuPbaa5GXl4e8vDykpaX59Py1S8Ff51D4MsZcelLHXBAEdHV1obGxUUzq\nXF9dl1tbWz0+bnh4+LCJnvMWRImJiYiKivL67zJ+zqUn55jL5sIGuvSqq6uxZcsWfPTRR6ioqBDX\nq1Qq5Obm4qabbsL8+fORkZHh9V90ROR/FAqFeDPkjIyMEbft6+tDc3PzsAmf62t3dze6u7tx4sSJ\nEY/pvJjEmdQ5i2vdYDBIepNrorFiT9w4YLVasWXLFrz99tsoLS0V10dERODWW2/FggULMHfuXERG\nRnqxlUREY+Mczh0uwXM+5aW+vt7j27NotVq3JG+oxC82NpZ/7NIlx+FUF+MpiausrMT69euxadMm\ndHZ2Ajh7r6gFCxYgPz8f8+bN47w2IhpXurq6UF9fLyZ1zuJab2pqgsPhGPVYQUFBg3rwXJO8pKQk\nJCYmslePLgiTOBdyTeLOZzy/tLQUzz77LHbu3CmumzVrFpYuXYpf/epXnADsITnPoZArxlx6jPlg\nDocDzc3NwyZ5zmWLxTLqsRQKBfR6PYxGI5KSkpCUlIS+vj7k5uaK63hz80tPzp9zzokbJw4cOIAn\nn3wSRUVFAM72ut1///1YtmwZsrOzvdw6IiJ5UCqV4n34ZsyYMex2VqvVbaj23FJbW4uGhgbxgo1D\nhw6J+65fv15cDggIQEJCglui5yzOdfHx8bzdCg3Cnjg/UF9fj9WrV+PDDz8EcHau28MPP4wVK1ZA\np9N5uXVEROOX3W5HY2Mj6urqUFtbKxbXelNT06jfSyqVCgaDYchELzk5GSkpKZyj56c4nOrCn5I4\nQRDw5ptv4vHHH4fVakVwcDD+9Kc/YdWqVXxKAhGRTPT19aGxsXHYJK+2thbNzc2jHkej0YgJnWtx\nrjMYDOzNkyEmcS7kmsSdO57f1taGBx98EJ999hkA4LbbbsOrr76KCRMmeKmF/kfOcyjkijGXHmMu\nvbHEvLe3VxyidU3yampqYDKZYDKZ0NXVNeIxVCoVjEbjsElecnKy3z4DV86fc86J8zPHjx/H9ddf\nj9raWkRFReHdd9/Fbbfd5u1mERHRJRIcHIy0tDSkpaUNu01nZ6eY0A1VzGYzqqurUV1dPewx9Hr9\nsEleSkoKR3l8DHviZOaHH37Addddh+bmZsycORNbtmxBSkqKt5tFREQ+rqenB7W1tcMmeXV1daM+\nAzciIgIpKSlITU1FamoqJkyY4FYiIiIkOpvxg8OpLuScxNXX1yMnJwctLS2YO3cutm7dCrVa7e1m\nERGRH3A4HGhsbByxN89ms414jJiYmEGJnbOkpKT47XDtpcQkzoVck7iioiI899xzKCoqwrXXXosd\nO3bwP8MlJuc5FHLFmEuPMZeeXGMuCALa2tpgMplQXV2Nqqoqt1JdXY2enp4Rj2EwGIZN8oxG4yW7\n8EKuMQdkNieusLAQK1asgMPhwEMPPYQ1a9YM2uaRRx7Brl27oFarsXHjRuTk5Hi8r1zt2rULRUVF\n0Ol0+PDDD5nAERGRpBQKBbRaLbRaLS6//PJB7wuCALPZPCi5c5aamhrxvnoHDhwYtL9KpUJycjLS\n0tIwceLEQSUsLEyK05QVn+qJczgcmDx5Mnbv3o3ExERcccUV2Lx5MzIzM8Vtdu7ciYKCAuzcuRMH\nDx7Eo48+irKyMo/2BeTZEycIAqZOnYrKykps2rQJCxcu9HaTiIiIzovdbkddXd2wSV5jY+OI+8fF\nxYkJ3bmJnl6v99t75F1wT9zs2bOxfPly3H333Zf0WZvl5eVIT09HamoqACA/Px/btm1zS8S2b9+O\nBx54AAAwc+ZM8aHHVVVVo+4rV4cOHUJlZSX0ej3uuusubzeHiIjovKlUKvGCiLy8vEHvnzlzBiaT\nCadOncJPP/3kVk6dOoXm5mY0NzejtLR00L4ajWZQYuesp6SkIDAwUIpTlJxHSVxwcDAWLVqElStX\n4v7778fSpUuRkZFx0RtTX1+PpKQksW40GnHw4MFRt3E+y260feXqhx9+AABcd911fJCyhOQ8h0Ku\nGHPpMebSY8yHFhoaioyMjCHzi4GBATQ0NAxK7JzL7e3t+P777/H9998P2lepVGLSpEkoKytDZGSk\nFKciGY+SuOLiYhw/fhxvvfUW/vGPf+CNN97A1VdfjWXLluGOO+64aBmup12hchsOvVAmkwkAxF5G\nIiKi8SQgIABGoxFGoxHXXHPNoPctFouY0B09ehR79uzB/v37AZydqnX8+HE8//zzWLdundRNv6Q8\nvrAhIyMDf/3rX/H888/jn//8JzZs2IB77rkHOp0OixYtwtKlS0e8CaEnEhMTUVtbK9Zra2thNBpH\n3Kaurg5GoxH9/f2j7uu0aNEiMSGKiorC9OnTxb+KiouLAcCn6i0tLQDOfkh9oT3jpZ6bm+tT7RkP\ndec6X2nPeKk7+Up7WGd9tHpHRwc+/vhj1NfXIyQkBKdOncKhQ4fQ0NCA1tbWITt7amtrUSyD3y/O\n5ZFuyuw05gsbDh8+jJUrV2Lfvn1nD6RQ4LbbbkNBQQHi4+PHckjY7XZMnjwZe/bsgcFgwJVXXjni\nhQ1lZWVYsWIFysrKPNrX2U659eR99tlnWLBgAX7+85+L8SYiIvJXQw2fupaOjo5h91UqlUhJSRk0\nNy43NxcxMTESnsXFcdFuMWKz2bB582asX78e33zzDSZPnozXX38dd955J3bs2IGnnnoK99xzD4qK\nisbUUJVKhYKCAtxwww1wOBxYvHgxMjMzsWHDBgDA0qVLceONN2Lnzp1IT0+HRqPBe++9N+K+/mDO\nnDkIDg7G/v37cezYMb85L1/n+hcbSYMxlx5jLj3G/CyLxeJ2jznn8k8//YSqqir09vYOu69Goxny\nKtWJEyciOTl50DSv4uJiWSZwo/Eoifvuu++wYcMGfPDBB7DZbLj11luxbt06XHvtteI2S5YsQXx8\nPO68884LatD8+fMxf/58t3VLly51qxcUFHi8rz+IiorC3Llz8fnnn+OJJ57AJ5984reXUhMRkX+w\nWq2DbgrsWu/s7Bxxf9dbipxb4uLi+D0ID4dTAwICYDAYsGTJEvz2t79FQkLCkNsdPXoUv/vd77B3\n796L3tCLRY7DqQBQXV2NadOmobu7G2+99RaWLFni7SYREdE41tvbO+TTG5x153zu4Wg0GkyYMMHt\nGaypqalIS0tDWloawsPDJToT33bBj9365JNP8Mtf/vKSPQ5DSnJN4gDggw8+wMKFC6FSqbBlyxbc\nfvvt3m4SERH5qXNvzuuarFVXV6OhoWHE79OgoCAxQXNN1Jx1rVbL3jQP8NmpLuSaxDnnUKxevRov\nv/wylEolNm7cyKc3XEKctyI9xlx6jLn0fCXmp0+fRk1NDUwmE2pqasTirNfV1cHhcAy7v1KpRFJS\nklti5pqoxcfHIyAgQMIzGp6vxHwsZPXsVBrZunXrEBgYiOeffx733Xcf9u/fj9deew2hoaHebhoR\nEfmIgYEBNDU1DUrMXOsjXeHpdO4D610TNaPRCJWKaYQ3sSdOhgRBwPr167Fy5Ur09vYiKysL69ev\nx1VXXeXtphERkQRsNptbUnZuolZbW4v+/v4RjxEcHIzk5GSkpKQgOTl50LLRaERISIhEZ0TD4XCq\nC39I4pwqKipw991343//+x8AYOHChVi3bh0MBoOXW0ZERGPlcDhgNptRV1c3bKLW2to66nF0Ot2Q\nyZlzWafTcU6aDDCJcyHXJG648XybzYYXX3wRL730Enp7e6FWq7F8+XKsWrVq2KuIyTNynkMhV4y5\n9Bhzadntdnz66adISkpCXV3dkKWhoQF2u33E4wQGBopJ2VCJWlJSEtRqtURn5fvk/DnnnDg/plar\n8cwzz2DRokVYtWoVtm7dildffRUFBQVYvHgxVq5cifT0dG83k4jI7/X19aGhoWHY5Kyurg6NjY0Y\nGBgY9Vg6nQ5GoxFJSUlD9qLp9XqfuWiAvIc9cX7m8OHDeO655/Dpp5+K6+bOnYtly5bhlltuGXQX\nayIiGp3VakVDQ8OISZrZbB71+0WhUECv14sJmvOh7q7FYDBwLhqJOJzqwt+TOKfKykq88sor+Oij\nj9DT0wMASEhIwL333ov8/Hz87Gc/41wIIhr3ent70djYKCZo55b6+no0NDSgq6tr1GMFBAQgISFh\n2OTMaDQiISEBQUFBEpwZ+QsmcS7kmsSNdTy/o6MDmzZtwvr163Hs2DFxfXp6OvLz83HXXXchKyuL\nCd0Q5DyHQq4Yc+n5a8z7+/thNpuHTc6cpa2tzaPjBQcHw2AwjJikxcfHe3TLDX+NuS+Tc8w5J24c\ni46OxiOPPILf//73KC0txebNm/Hxxx/j5MmTePbZZ/Hss88iJSUFN998M26++Wbk5uayG5+IfFZ/\nfz+am5tHTdCam5s9+oNdpVIhISEBBoNhxBIdHc0/dsnnsCduHLLb7SgpKcHmzZuxfft2t+fbqdVq\n5OXl4dprr0VeXh6ys7M5eZaILimHw4HW1laYzWY0NTWhqalp2GVPe86cc89GS850Oh1/x5FP43Cq\nCyZx7gYGBnDo0CF8/vnn2LFjBw4fPuz2fnR0NK655hrk5eXhqquuwrRp03hxBBGNShAEtLe3j5qY\nmc1mNDc3e3TFJnB23plOpxs1QdPr9XyaAPkFJnEu5JrESTWe39DQgD179qCoqAh79+6FyWRyez80\nNBQzZszA7NmzMWvWLMyaNctv70cn5zkUcsWYS+98Yj4wMID29nY0NzejpaUFZrPZLRlzTdDMZvOo\nTwxwFRsbi/j4eOj1esTHx7stu67TarVQKpVjPFvfwM+59OQcc1nMiWtvb8fdd98Nk8mE1NRUfPzx\nx4iKihq0XWFhIVasWAGHw4GHHnoIa9asAQA8/fTTePvtt6HT6QAAL7zwAubNmyfpOfgDg8GA++67\nD/fddx8AoKqqCnv37kVJSQlKS0tx4sQJ7Nu3D/v27RP3SUhIwPTp05GTk4OcnBxMnz4daWlpHKIg\n8nEDAwPo7OzE0aNH0dLSgpaWFjFBG2q5ra3N4x4zAIiMjPQoMdPpdLxik2gMfKYnbvXq1dBqtVi9\nejXWrVuHjo4OvPjii27bOBwOTJ48Gbt370ZiYiKuuOIKbN68GZmZmVi7di3Cw8Pxhz/8YcSfI9ee\nOF/R2tqKgwcPoqysDGVlZSgvLx/y0vvw8HBMnToVU6ZMcStGo5GTg4kuEZvNhtbWVrS1tYmltbUV\nra2tFyUpA85OsdDpdIiLi4NOp0NCQsKwSRovkiK6cLIYTs3IyEBJSQn0ej2ampqQm5uL48ePu21T\nWlqKtWvXorCwEADEJO+xxx7D2rVrERYWhlWrVo34c5jEXVwDAwOoqqrCkSNHUFFRgSNHjuDIkSNo\nbGwccvvw8HBkZmbisssuQ3p6uluJiYlhgkeEs/PJOjs7h0zIRqo77wl5Ps5NykZa1mq1nBNLJDFZ\nDKeazWbo9XoAgF6vh9lsHrRNfX09kpKSxLrRaMTBgwfF+ptvvon3338fM2bMwKuvvjrkcKxc+ep4\nfkBAACZOnIiJEyfizjvvFNebzWYcPXoUlZWVOHr0qFhaWlpQXl6O8vLyQceKiopCeno6JkyYgJSU\nFLeSnJws+b+nr8bcn/lTzAVBgM1mQ0dHh1gsFotb/dx1zmSsvb0dDofjvH9mSEgIYmNjxaLVasXX\n4ZKyAwcO+E3M5cKfPudy4a8xlzSJmzt3Lpqamgatf+6559zqCoViyB6ZkXppli9fjieffBIA8H//\n939YtWoV3nnnnSG3XbRoEVJTUwGcTRymT58u/uMWFxcDgM/VnXylPZ7U9Xo9FAoFsrKyxPe3bt0K\nk8mE8PBwnDx5El999RXq6+vR1NQEi8WCQ4cO4dChQxhKRESE+IU0depUGAwGdHd3Q6vV4vrrr4fB\nYMDx48ehUql84vxZP/96RUWFz7RHEAQUFhbCarUiKysLXV1d2LdvH6xWKxITE9HR0YGKigp0d3dD\no9Ggo6MDJpMJ3d3d6O3thcViOa+J/eeKiIiAWq1GREQE0tLSEBsbizNnziAiIgIzZsyAVqtFXV0d\nIiIicMMNN0Cr1Yp/HI12ftOmTRPrFRUVPhHv8VR38pX2sO5bdedydXU1RuNTw6nFxcWIj49HY2Mj\n8vLyBg2nlpWV4emnnxaHU1944QUEBASIFzc4VVdX45ZbbsH3338/6OdwONU3CYKAlpYWnDhxAtXV\n1TCZTGKpqamByWSCzWYb9TgKhQIxMTFuPQ7nvmq1WkRHRyMmJgbR0dEIDw/nMK6fEAQBfX19sFqt\nOH36NLq6usZcuru7z3u+2LlCQ0MRHR2NqKgoREdHu5Wh1sXExCA2NhYxMTGc6E9EAGQyJ2716tWI\njY3FmjVr8OKLL8JisQy6sMFut2Py5MnYs2cPDAYDrrzySvHChsbGRvFWF6+99hq+/vprfPjhh4N+\nDpM4eRIEAW1tbaipqXG7K/u5zzz05AHU51IqlW5fotHR0YiMjERERATCw8MREREx5HJYWBjUajXU\najU0Gg3UajUCAwOZEA5BEATY7Xb09fWht7fXrdhsNthsNlitVrG41s/3vbEMQw7H2RPmWsLDwwcl\nX8MlZsHBwRetLUQ0PskiiWtvb8ddd92Fmpoat1uMNDQ0YMmSJdixYwcAYNeuXeItRhYvXow///nP\nAID7778fFRUVUCgUmDBhAjZs2CDOsXMl1ySu2E/H8y82u93udh8r1yvynK/OOUfOeUinT5++aD9f\nqVSKiZ1arUZoaCiCg4NHLYGBgQgICIBSqRRfh1se7X3nPbQGBgbgcDhGffVkm4GBAfT394uJ11DJ\n2GhFqv93KpUKGo3GLekeSwkPDx8XN4vl7xbpMebSk3PMZXFhQ0xMDHbv3j1ovcFgEBM4AJg/fz7m\nz58/aLv333//kraP5EGlUiEuLg5xcXEe79PX1weLxSImdu3t7ejq6sLXX38Ng8HgNrzmfO3s7BR7\ngZw9STabDf39/eju7kZ3d/clPEt5CggIGDKBDQ0NhUajgUajgc1mQ0pKitiz6Vx/bn24ZY1Gw6sn\niWjc8JmeOKnItSeO5KG/v98tqbPZbB71XNnt9mF7xUbrNRtqWRCEIXvvRnsd7T1PehVdS1BQkLg8\nHnq1iIguNlkMp0qFSRwRERHJxUh5S4DEbaExOvfSdLr0GHPpMebSY8ylx5hLz19jziSOiIiISIY4\nnEpERETkozicSkRERORnmMTJhL+O5/syxlx6jLn0GHPpMebS89eYM4mTCeczJUk6jLn0GHPpMebS\nY8yl568xZxInExaLxdtNGHcYc+kx5tJjzKXHmEvPX2POJI6IiIhIhpjEyUR1dbW3mzDuMObSY8yl\nx5hLjzGXnr/GfNzdYiQ3NxclJSXebgYRERHRqK655pphL8wYd0kcERERkT/gcCoRERGRDDGJIyIi\nIpIhJnE+prCwEBkZGZg0aRLWrVs35DaPPPIIJk2ahOzsbBw5ckTiFvqf0WL+wQcfIDs7G9OmTcNV\nV12F7777zgut9C+efM4B4Ouvv4ZKpcKnn34qYev8kycxLy4uRk5ODrKyspCbmyttA/3QaDFvbW3F\nvHnzMH36dGRlZWHjxo3SN9KP/OY3v4Fer8fUqVOH3cbvvj8F8hl2u12YOHGiUFVVJfT19QnZ2dnC\n0aNH3bbZsWOHMH/+fEEQBKGsrEyYOXOmN5rqNzyJ+VdffSVYLBZBEARh165djPkF8iTmzu3y8vKE\nm266SfjXv/7lhZb6D09i3tHRIUyZMkWora0VBEEQWlpavNFUv+FJzJ966inhscceEwThbLxjYmKE\n/v5+bzSjh141AAAFMElEQVTXL/z3v/8VDh8+LGRlZQ35vj9+f7InzoeUl5cjPT0dqampCAwMRH5+\nPrZt2+a2zfbt2/HAAw8AAGbOnAmLxQKz2eyN5voFT2I+e/ZsREZGAjgb87q6Om801W94EnMAePPN\nN3HnnXdCp9N5oZX+xZOYf/jhh7jjjjtgNBoBAFqt1htN9RuexDwhIQFdXV0AgK6uLsTGxkKlUnmj\nuX5hzpw5iI6OHvZ9f/z+ZBLnQ+rr65GUlCTWjUYj6uvrR92GScXYeRJzV++88w5uvPFGKZrmtzz9\nnG/btg3Lly8HACgUCknb6G88ifmJEyfQ3t6OvLw8zJgxA5s2bZK6mX7Fk5gvWbIElZWVMBgMyM7O\nxhtvvCF1M8cVf/z+ZMrvQzz9ohLOuSsMv+DG7nxit3fvXrz77rs4cODAJWyR//Mk5itWrMCLL74I\nhUIBQRAGfebp/HgS8/7+fhw+fBh79uyBzWbD7NmzMWvWLEyaNEmCFvofT2L+/PPPY/r06SguLsZP\nP/2EuXPn4ttvv0V4eLgELRyf/O37k0mcD0lMTERtba1Yr62tFYc2htumrq4OiYmJkrXR33gScwD4\n7rvvsGTJEhQWFo7YXU+j8yTm33zzDfLz8wGcnfy9a9cuBAYGYsGCBZK21V94EvOkpCRotVqEhoYi\nNDQUV199Nb799lsmcWPkScy/+uorPPHEEwCAiRMnYsKECfjxxx8xY8YMSds6Xvjl96d3p+SRq/7+\nfiEtLU2oqqoSent7R72wobS01C8mZnqTJzE3mUzCxIkThdLSUi+10r94EnNXixYtEj755BMJW+h/\nPIn5sWPHhOuuu06w2+2C1WoVsrKyhMrKSi+1WP48ifnKlSuFp59+WhAEQWhqahISExOFtrY2bzTX\nb1RVVXl0YYO/fH+yJ86HqFQqFBQU4IYbboDD4cDixYuRmZmJDRs2AACWLl2KG2+8ETt37kR6ejo0\nGg3ee+89L7da3jyJ+TPPPIOOjg5xflZgYCDKy8u92WxZ8yTmdHF5EvOMjAzMmzcP06ZNQ0BAAJYs\nWYIpU6Z4ueXy5UnMH3/8cTz44IPIzs7GwMAAXnrpJcTExHi55fL161//GiUlJWhtbUVSUhLWrl2L\n/v5+AP77/cnHbhERERHJEK9OJSIiIpIhJnFEREREMsQkjoiIiEiGmMQRERERyRCTOCIiIiIZYhJH\nREREJENM4oiIiIhkiEkcERERkQwxiSMiIiKSISZxRERjYLVakZGRgZkzZ8Jut4vrv/zySwQEBODv\nf/+7F1tHROMBH7tFRDRGFRUVmDVrFlauXIkXXngBZrMZ2dnZmD17Nv797397u3lE5OeYxBERXYDX\nX38df/zjH/HFF1/g5ZdfRmVlJb799ls+yJyILjkmcUREF+imm27Cnj17YLfb8Z///Ad5eXnebhIR\njQOcE0dEdIEWLlyIvr4+ZGdnM4EjIskwiSMiugBNTU149NFHcfnll6OiogJ/+9vfvN0kIhonmMQR\nEY2RIAh44IEHEBoait27d2PFihVYs2YNvv/+e283jYjGAc6JIyIao1deeQWPPfYY9u7dizlz5qC/\nvx+zZs1Cb28vDh06hJCQEG83kYj8GHviiIjG4PDhw/jLX/6Cxx9/HHPmzAEABAYGYvPmzaiursaq\nVau83EIi8nfsiSMiIiKSIfbEEREREckQkzgiIiIiGWISR0RERCRDTOKIiIiIZIhJHBEREZEMMYkj\nIiIikiEmcUREREQyxCSOiIiISIaYxBERERHJ0P8DTXjYHxEkAyYAAAAASUVORK5CYII=\n", "text": [ "" ] } ], "prompt_number": 2 }, { "cell_type": "markdown", "metadata": {}, "source": [ "The contour defining the airfoil will be partitioned into `N` panels, using the same method as in the [Lesson 10](http://nbviewer.ipython.org/github/barbagroup/AeroPython/blob/master/lessons/10_Lesson10_sourcePanelMethod.ipynb)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We define a class `Panel` that will store all information about one panel: start and end points, center point, length, orientation, source strength, tangential velocity and pressure coefficient. We don't save the vortex-sheet strength because all panels will have the same value." ] }, { "cell_type": "code", "collapsed": false, "input": [ "class Panel:\n", " \"\"\"Contains information related to one panel.\"\"\"\n", " def __init__(self, xa, ya, xb, yb):\n", " \"\"\"Creates a panel.\n", " \n", " Arguments\n", " ---------\n", " xa, ya -- Cartesian coordinates of the first end-point.\n", " xb, yb -- Cartesian coordinates of the second end-point.\n", " \"\"\"\n", " self.xa, self.ya = xa, ya\n", " self.xb, self.yb = xb, yb\n", " \n", " self.xc, self.yc = (xa+xb)/2, (ya+yb)/2 # control-point (center-point)\n", " self.length = sqrt((xb-xa)**2+(yb-ya)**2) # length of the panel\n", " \n", " # orientation of the panel (angle between x-axis and panel's normal)\n", " if xb-xa <= 0.:\n", " self.beta = acos((yb-ya)/self.length)\n", " elif xb-xa > 0.:\n", " self.beta = pi + acos(-(yb-ya)/self.length)\n", " \n", " # location of the panel\n", " if self.beta <= pi:\n", " self.loc = 'extrados'\n", " else:\n", " self.loc = 'intrados'\n", " \n", " self.sigma = 0. # source strength\n", " self.vt = 0. # tangential velocity\n", " self.cp = 0. # pressure coefficient" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 3 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Like before, we call the function `define_panels()` to discretize the airfoil geometry in `N` panels. The function will return a NumPy array of `N` objects of the type `Panel`." ] }, { "cell_type": "code", "collapsed": false, "input": [ "def define_panels(x, y, N=40):\n", " \"\"\"Discretizes the geometry into panels using 'cosine' method.\n", " \n", " Arguments\n", " ---------\n", " x, y -- Cartesian coordinates of the geometry (1D arrays).\n", " N - number of panels (default 40).\n", " \n", " Returns\n", " -------\n", " panels -- Numpy array of panels.\n", " \"\"\"\n", " R = (x.max()-x.min())/2 # radius of the circle\n", " x_center = (x.max()+x.min())/2 # x-coord of the center\n", " x_circle = x_center + R*np.cos(np.linspace(0, 2*pi, N+1)) # x-coord of the circle points\n", " \n", " x_ends = np.copy(x_circle) # projection of the x-coord on the surface\n", " y_ends = np.empty_like(x_ends) # initialization of the y-coord Numpy array\n", "\n", " x, y = np.append(x, x[0]), np.append(y, y[0]) # extend arrays using np.append\n", " \n", " # computes the y-coordinate of end-points\n", " I = 0\n", " for i in xrange(N):\n", " while I < len(x)-1:\n", " if (x[I] <= x_ends[i] <= x[I+1]) or (x[I+1] <= x_ends[i] <= x[I]):\n", " break\n", " else:\n", " I += 1\n", " a = (y[I+1]-y[I])/(x[I+1]-x[I])\n", " b = y[I+1] - a*x[I+1]\n", " y_ends[i] = a*x_ends[i] + b\n", " y_ends[N] = y_ends[0]\n", " \n", " panels = np.empty(N, dtype=object)\n", " for i in xrange(N):\n", " panels[i] = Panel(x_ends[i], y_ends[i], x_ends[i+1], y_ends[i+1])\n", " \n", " return panels" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 4 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we can use our new function to define the geometry for the airfoil panels, and then plot the panel nodes on the geometry." ] }, { "cell_type": "code", "collapsed": false, "input": [ "N = 40 # number of panels\n", "panels = define_panels(x, y, N) # discretizes of the geometry into panels\n", "\n", "# plots the geometry and the panels\n", "val_x, val_y = 0.1, 0.2\n", "x_min, x_max = min( panel.xa for panel in panels ), max( panel.xa for panel in panels )\n", "y_min, y_max = min( panel.ya for panel in panels ), max( panel.ya for panel in panels )\n", "x_start, x_end = x_min-val_x*(x_max-x_min), x_max+val_x*(x_max-x_min)\n", "y_start, y_end = y_min-val_y*(y_max-y_min), y_max+val_y*(y_max-y_min)\n", "\n", "size = 10\n", "plt.figure(figsize=(size, (y_end-y_start)/(x_end-x_start)*size))\n", "plt.grid(True)\n", "plt.xlabel('x', fontsize=16)\n", "plt.ylabel('y', fontsize=16)\n", "plt.xlim(x_start, x_end)\n", "plt.ylim(y_start, y_end)\n", "plt.plot(x, y, color='k', linestyle='-', linewidth=2)\n", "plt.plot(np.append([panel.xa for panel in panels], panels[0].xa), \n", " np.append([panel.ya for panel in panels], panels[0].ya), \n", " linestyle='-', linewidth=1, marker='o', markersize=6, color='#CD2305');" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAnEAAAB+CAYAAABLREfEAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xt8U/X9+PHXyb33e5M2KS20XEVBxQv6RUG8DSeb3zl1\nuq+XCUPnDXQKsk3Efad4mVc2ZbLp/Kns7hcnl00Z1Tm5iIo3vBRpS3pL27TpJW2TnOT8/igNlAQo\nSNMmvJ+Px3kkJzlJ3nn3NHnn8/mcz1E0TdMQQgghhBBxRTfUAQghhBBCiMMnRZwQQgghRBySIk4I\nIYQQIg5JESeEEEIIEYekiBNCCCGEiENSxAkhhBBCxCHDUAcQa9OnT+fNN98c6jCEEEIIIQ7p7LPP\npry8POp9yrE2T5yiKMTjW7722mt5/vnnhzqMY4rkPPYk57EnOY89yXnsxXPOD1a3SHeqEEIIIUQc\nkiIuTpSUlAx1CMccyXnsSc5jT3Iee5Lz2EvUnEsRFyemT58+1CEccyTnsSc5jz3JeexJzmMvUXMu\nRZwQQgghRBySIk4IIYQQIg7J0alCCCGEEMPUweqWY26eOCHixca1a3j1mafQ+X2ETGZm33ALM2Zd\ndMzGIYQQoj8p4uJEeXl5wg7MHK7Ky8vRurxDUsBsXLuGPy9ewJzOqvBtK+/+Ck3TjvrrK4pyeHEs\n3gUwKHkYyH4uReXRJZ8tsSc5j71EzbkUcUIcwAebN1Hxh+f6FTDPLKygvq6OcSedjNfrpaurK3zZ\nt+y73t3dTU9PD36fj0BPD6GeboK+bkI9PWj+HvD50Pw+lGAAfVBFHwyiD6m4Ozp5OK9/cTXHW83C\ny2fzvEGhr2Fd22cJr2v9b9t/u75m+b51nU6HTq9Hp9Oj0+vR63svdXo9dZ52HswJ9Y+js4pF113J\ni+MmoBmMYDKByYJiMoPJhM5kwWixkJScTFJSUsRisVii3p6UlITb7cbj8WCxWDCbzREFZqyLSiGE\nGM6G3Zi49evXM3/+fILBIHPmzGHhwoUR29x6662sW7eO5ORknn/+eU488USgdx6Y9PR09Ho9RqOR\nrVu3RjxWxsTFvyNtiVFVlebmZtxuNy0tLbS2ttLa2hq+3tLSQru7mZ7mRrS2VtxffckjWYGI57nH\npXJikg6TAuY9i2nfS50Svr7v7Rrg03oXf2jPpQY+Tdvneu/ltq4Qi/Mjf2MtbVQpNe0tbJT9FpR9\nrhPluqL0uy3iOZT+t23vjh7HI00qF6bpsChK7/vU9c+DAhHvqfd9a3tzsP/9of55CAAhvRHNaAKj\nCUxmGls9PBTlb7LEWMQZV88hLS2N1NTUg15GKw6/DmkZFEIMprgZExcMBrn55pt54403sNvtnHLK\nKcyePZvx48eHt1m7di07d+6koqKCLVu2cOONN7J582ag942Wl5eTnZ09VG9BDLJoLTFP3/UFn33+\nOdbiEurr62lsbKSxsRF3owuvqx5fcyMhjxudt4NMvUKmDjL0hK/n6KFUr5CpB4sCniB4QvB/apBo\n/yI6g4G2vAL0liR0Zgv6pGQMScnok5IxJqdgSk3FlJKGOTUNU2oaxuQUjMkpmJOTMZvNmM1mUsxm\nsvdcN5vNmEwmDAZDeKm57kqoeDvitUtnXMCjq9cdtXxqmoaqqqiqSiAQ6LeoqsruuVdD1ZaIx5mO\nn0L+T+8LP64zEKB1z2MCgQA9XV78nR34vV4C3k4CXZ0EurwEu7pQe7oIdncT8vWg+XrA39saScCP\nEvBjUAOYVBWDFsSsBDAFVUyhLsx+eFcNRf2bFDRUMu3Ze/CG2LNoNPVd13rX++7rQiFosqBZksCS\ngpKSij41HXNmNpasbFKzc8jIzCQjI4PMA1xmZGRgMBikZVAIMaSGVRG3detWysrKwjMrX3HFFaxe\nvbpfEffqq69yzTXXAHDaaafh8XhwuVxYrVaAhG1lS9T+/IPx+/04nU6qq6vZvXs31dXVvPP73/CA\nrqnfdjf21HLvkjuZmaojUw9FeoUMHSTroC0IbSHwGDTa0vR0G80ELMkEk9Pwp2XQnJFNW3Yurnwb\nydYCMmyFZGVnU5iRQdvNN0DnlxFxlZ0186gWUtF855Y7WLm4pl9x8GxqCZfNu/movo6iKBiNRoxG\nI0lJSRH3f+/2Razcr0h5NrWEaxYvYcaFFx7VWKD/fh4MBunu7u63fDbvOqjZFvE4l72MdedfhL+t\nFX+bh2BnG6HODkJdndDdha6nC73fR2bAh00LkkIPKf4eUtRWUroUUpohRde7GBT2KQahNaRRo9Gv\nQPSGIGAw4fQFeTCvfyxzOqu45/Zb2PzhR2RmZpKdnR2xpKenH9XWwK/jWPxsGWqS89hL1JwPqyKu\ntraWoqKi8LrD4WDLli2H3Ka2thar1YqiKJx77rno9XrmzZvH3LlzYxa72Gug3UuaptHU1MTOnTv3\nLhUVNO7aiVpbRVJbCw4jFBkVHAaYalSoaQtCTuRuq5gt1I6fRIu1gNQCB2n2InIcI8izWpmQn09e\nXh7Z2dno9foBv49zr7iKlfuNiRuMQiqavny9uGI5iq8HzWzhsnk3x7x1Zyjj0Ov1pKamkpqaGr7t\nqrt+ErWonH//YwOOye/34/V66ejooLOzk46ODtrb22lra2N3WxttbjdedxPd7iZ8nhZ8nlbUjnaC\nne3Q7UXp6ULv6yEz5KdN04j2MVroquKM39xDWwjaghpfhnp/ULSHNNqD0IGCaklBS01Hl56JMTMH\nS24+yfk2snJzoxZ+2dnZZGRkHHIflu5dIY4dw6qIG+gv0wO1tr399tsUFhbS1NTEeeedx7hx45g2\nbdrRDHHIxMsviGjdS8/evZOq6mpScnLZsWMHn+3YgfvLHfidleQEunEYFRxGhXFGhXMNEAJq0KhJ\n1dFuSaUjx8qXhSOoKx2Hp3wDdFREvG7pmWfz0FFuHVuw6G42nnDCkBVSM2ZdNCy+fGMZx6H286NR\nVJpMJkwmE1lZWUccp6ZpdHZ2ctd3LoZd70Tc32gv4/XzvoGvuRFfSzOqx02o3YPS2YG+24tV9ZEe\n9JLR4SW9q56MRoWMryBNBz2h3tbj9hDsCmps31MItoWgIwQBczJaahr69CxM2bkk5dtIsxWSnW+l\nraGO2jV/Yb7WHI5l5d1f9cvd/uLlsyWRSM5jL1FzPqyKOLvdjtPpDK87nU4cDsdBt6mpqcFutwNQ\nWFgIQF5eHpdccglbt26NWsRde+214S7bzMxMJk+eHP4Dl5eXA8j6EaxrmsZjP1vMN5t2QVLvyUC2\ndYc4sbuKv9x1M+en6lCCGtP0UGbWUWPR2Kwz0J6aSc/oMdSXjeNTSyr5xSP59re/jd1u5+233+73\neo8te4ClK3/FEhrDz/9ako0Fe1rHjvb7U5JTmL1gYb/7922WH075P1bWleSUcHd23/19YhlPWloa\nY865gKW7dkbujw/1tgwe6PFnnHEGHo+HdevW0d7eTm5JCTUtLWzdshlvaytWixGfu4nqykqCXZ3Y\nCaDv9tLj7SbV18UkpZv07kacNZ+TooOzUnT4NLivUeW7Gfp+/3+Tuyt54Hvf5vmCAlw6E8aMLEZN\nmEiarRBPWxsZGRmcddZZ5ObmsnPnTjIzM7ngggtink9Zl3VZ713vu15VVcWhDKujU1VVZezYsWzY\nsIHCwkJOPfVUVq1aFXFgw/Lly1m7di2bN29m/vz5bN68ma6uLoLBIGlpaXi9Xs4//3yWLFnC+eef\n3+814vXo1PJh1p8fCoXCB5d88MEH7Ny2Be2Lj9B3d7HEGuWoyhaFsZNOJGv88ThOPpUxE0+grKyM\n7Ozswx4btHHtGv6+T0vMxYPUOjbccn4siMecx2p/hN5xgh6PB7fbjdvtprm5mebmZpoaG/E01PHp\n31axxOiJeNz/NqqckawjSw9ZeoVsfe+YUU8Q3vaGKDAqtAahNajREoQuoxktLQN9dh6mXCspBXay\nCh3kW63k5+dj3XOZn59PdnY2Op3uoPmR7t3+4nE/j3fxnPO4OTrVYDCwfPlyLrjgAoLBINdffz3j\nx49nxYoVAMybN49Zs2axdu1aysrKSElJ4bnnngOgoaGB//7v/wZ6i8GrrroqooATR66trY0tW7aw\nefNmtrzzH5re20yJv5MTLAozLAqzdfCxXmOdLvp4ndJp5xy17s7h0s0oBMR2f9Tr9eTk5JCTkxP1\n/gVffQE7/hVxe+GZ07ngoSfCRd+OpibcLhed9TXs/PQTmhWVkKcFfWcbqV2d2PGT1dVElq+J7MbP\nyPyid/qYlmBv4ffJnmLPE9TwhBQCKWko6Vnoc/JIyi8gpbCInEI77a4Gal/7M7dpew9GkqN3hTh6\nhlVLXCzEa0vcYNv/1/K5V1+PZkli48aNbN3wOlR8yvEWmGRRGG9WaFShQmchUDKW9FP+izFnn8vk\nk07i8w/e5y8/uT3yYIDDGHguhDgyUcekHub/n6ZptLe309jYiMvlCk/Z01RbQ0etk25XHf6mBkKe\nFnSdbVh83WTrlT2tfISva8DTLUFuz41sK1jYaqB4wkQsNjtpI0aSZ3dQUFCAzWbDZrNRUFBARkbG\nsDmCV4ihdLC6RYo4wca1a/jT3fOZ660O3/bzRpUcncKsdB25BvjMr9CcbcN8/CmUnHcRp55zLkVF\nRVE/ZGPZvSSE6C/W/39+v5+mpqb+RZ/Lhbuuhk//uop7TZHdu8uaVGal6cjVK+ToIaBBcxCagxrN\nau9lm2IgmJaJkpWLyVpAcuEIskeUYCso6FfwWa1WjEbjIeOUbl0Rr6SI20e8FnGD0Z/f1tbG+vXr\nWbloAQ/qmyPuv6fNyNmXXcVJ376UM/5rWtR5xBJZPI+hiFeS89gbzJwvmH0h/xOle/fpEafynbt+\nSkNDA/V1dbTU7KbTWYW/sZ5QSxP69lbSgv5wkZdrUMjVg0UHLXuKvOYguNXey56kFMjIQZ9rJbnQ\nQVpRCVa7I1zo1VR8wdtPPcwPu3aHY1iZWsJ3h6iHQPbz2IvnnMfNmDgx+LxeL6+99hp/WLUK54a1\nnG8JMkLVos69NvbkU7nzV78ZgiiFEIlg9g23sHLxroju3St/fPchJ4vu7OykoaGht9Crr+fjhgZc\nNU46nFV019WgNrvA48bc3cmIQBc5rV3ktjvJ2f0e2e/2TszcHIQdqsY/OkMs2e/0cXM6q1hy+y3s\nqKzCbrdjt9spLCzEarViMMhXo4gP0hJ3DNA0jY0bN7Jy5UreXbOaGfoeLkrTEdTgk5wRfOiHn/ic\nEY978bjBPzOBECKxDXb3rqqqNDU19Sv4GurraK2uwlu3G7+rjvqdFfw8J/Jz/2m3yhWZehpVaFI1\nGoPQFAR/Sjpk52O0FZJaNJLcklHYHQ4KCwvDBd/hnHVDunLF1yHdqftI9CJu3w8Ln6JD5xjJm+Xl\njHLt4qI0HSNNCh+m5JMx+3tceNMC7Hb7URkMLYQQw9WBunWXmosZecrpeGuqUV11KK3NWLo7yNcr\n5Bsgz6CQrwej0tuq16hqvQVfUMOjMxHMyEafZyOpcAQZI0uxOYrCRV5fy96///mPyPPrDmFXrog/\nUsTtI16LuIH050crxh5u6p0fSrEkYTjnYs5fuIRRY8ZGfawcjNBfPI+hiFeS89g7FnJ+OD9U/X4/\n9fX11NXVUVtbS21tLQ3VVXRU76KnbjehZhf6thayQoHeIm9PsZej7z2jRrhVT9VoDinsCup4MLd/\nPNu6Q2wuO4Wfv/wXbDbbYZ0OUByZeN7PZUzcMeKPjz3ID/f5kAK4M8/AY9njee6d9w46zkPmXhNC\nJKrDOV2byWSiuLiY4uLiAz5f3zQsfUVeXV0d251OWqq+wuusJNBQB61NWLwdpIdUon3Vdm3fwjMn\nFtMUUvClZqLk2bDYR5A+cjSFI0dRVFREUVERDoeD/Pz8g06oLI5d0hKXALxeL/fffz//fmIZj9si\nf9G9VHYmv1y/cQgiE0KIY5ff7+fWi87jh5WbIu5b6DHhsBhI7u7AalCwGcC6p/u2W4MGFVyqhkvV\naNb0BNKz0eUXYHEUkzlqNI7iknChV1xcTFZWFoqiyPi7BCQtcQls9erV3HbzTZznrafEEP2PrJkt\nMY5KCCGEyWTi8gV3sTJKV+7iJ3u7cnt6eqitraWmpgan08mO3btp+qoC7+5dBOpr0LU0kurvwupt\nxupsxlb/CTkf9J45w6VqfK5Cg6rh0ZtoNqVi6vGyOEMNv9ZTt+/A7f5fLrnyKum2TUDSEhcn9u/P\n7+np4c4772TNiuUszddDajod37uJT175oxygcJTE8xiKeCU5jz3J+eDbf8xx0dSzWLDo7gE/vqur\ni9raWpxOJ06nk5rqalp27VvoNZEZ9LHTr/Gz/Mi2mXtdKg6LAV9aFjprIckjRpE1ejwjRo2iuLiY\nESNGMGLECMxm84DeSzy29MXzfi4tcQmi759H7exg28efkNPTwbN2I03nXMKlK17CYDCw8ZTTBzTu\nQwghRGzsP+a4vLz8sB6fnJzM6NGjGT16dNT7NU2jra2NhRefB/UfRtyv6vVowSDFXjcFu90U1H1C\n7rZX8QRht6qxJQD1qoY3OR3ybFgcJWSUjcNRWhYu8oqLi3n/P29HHmkr58IdUtISFyeiHV31mDvE\ntPl388MlPx+6wIQQQgwLB5pK5cXjZvKLP/wNp9NJdXU11dXV7K6sxL3zc7qqdxF01WJub8GqhwID\nFOwZo9ejQb0K9QGNelXjY7/CsvzIAyyeHnEKj7yyloyMjFi8zWOOtMQlgFefeapfAQewIEfHi+9v\nG5qAhBBCDCsHOkPGZfNuJikpiTFjxjBmzJioj1VVlbq6ut4Cb/duyquqaKz4gs6qnQTqnRhamkjB\nD0QWcW3vbWJpaTYeYxJaXgGWopFkjJnAiNFjGDlyJCNHjqSkpASLpf/47Hjtmh1OpIiLE7UNDVFv\nV3w9MY7k2BHPYyjileQ89iTnsTdYOT+cqVT2ZzAYwmPjotE0jZtnzYSKtyPuazWYSVM0rGo3hc2V\nFHqqKPhsI10a1AU0tu/pru1MSUfJt5NUXEqP0UTP1n9zu64l/DyD2TWbqPv5sCvi1q9fz/z58wkG\ng8yZM4eFCxdGbHPrrbeybt06kpOTef755znxxBMH/Nh45e7ojHq7HHkqhBCiz2DN+akoCpfecgcr\nF9dEPdJ2+jdm0djYSGVlJZWVlZTv+grXF5/RuetL1Hon5jY3jp4OCus+p7DpC/6vPcjtOZHns130\ng6t45fxvUjj+OErLyigtLaW0tLRfV6204O01oDFxU6dO5cYbb+Tyyy8f0NErRyoYDDJ27FjeeOMN\n7HY7p5xyCqtWrWL8+PHhbdauXcvy5ctZu3YtW7Zs4bbbbmPz5s0DeizE75i4kyaMp7T2C+7O27vT\ny5GnQgghYulIz+6jqiq1tbXhIu+1ZUv5SbAuYrtHm1Wuy9JjUaBWhZqARm1Aw2NKAWshHclpKLt3\nsdDcEX7MypRivvvA4/3iSKRC72ufdmv69Om89dZbZGVlcfXVVzNv3jzGjRt31APdtGkTS5cuZf36\n9QAsW7YMgEWLFoW3ueGGG5gxYwaXX345AOPGjaO8vJzKyspDPhbir4jbuHYNqx55gO2b/oPDqGAq\nm4A9P09OjSWEECJuHeggjMetkzj9quuo/vwzPF98is+5C31TA/kEsBsVNnlD3JEX2Ym4qFXPqJNO\nJWPMBAKahvv1v3Ob1sT73SE2d2nUBhX8SSlcfNMC5t9zbwze4dFzsLplQOfxKC8vZ8eOHVxzzTW8\n8MILTJgwgenTp/OHP/yBQCBw1AKtra2lqKgovO5wOKitrR3QNnV1dYd8bLzpOyL1ht1bmZOt56f5\nBrJDfr556495dPU6KeAG2eFOAyC+Psl57EnOY09yvucgjNSSfrc9m1rCdT+5lx/96Ec8+ORTrPjH\nv3h+RxUrG7u59aMaJv+pHEZFnvsbIFP1cdbn73DRupW4X/4N07pc/Nylsr4jxI9y9FySBg+ndfP2\nEw/w+H33Dvr7i5UBn4xt3LhxPProo9TW1vL73/8eVVW58sorcTgcLFy4kF27dn3tYBRFGdB28dSS\n9nVEOyJ1TmcVf1+xfGgCEkIIIY6CGbMu4rv3P8aLx83kpbIzefG4mQccHqQoCgUFBfzXf/0X+UXF\nvN8d4tfuICtagvzaHeT97hAZJ0/F/Pgf2Dzvf3EmZbCpSyPHoLB4v8mPF2dp/PPZX8XqbQ66wz6w\nwWKx8D//8z8cd9xxLFiwgH//+988/PDDPPLII1xyySUsX74cm812RMHY7XacTmd43el04nA4DrpN\nTU0NDoeDQCBwyMf2ufbaaykpKQEgMzOTyZMnh49a6fuFNBzWdX4f27pDAExJ6q23t3WHqKmvD7+X\n4RRvoq1Pnz59WMVzLKz33TZc4jlW1vsMl3hk/dhYV5JTmL1gYb/7D/X/H8ixsb5Hz+IcLfz9uLZL\nAZ2JF3+9nO46J4H2dm4q0LPEpbKtO8SUJB1TknTh7U2aOize/8H+H8vLy6mqquJQDmuy366uLlat\nWsUzzzzDe++9x9ixY7nxxhu59NJLWbNmDUuWLGHcuHH861+R/dwDoaoqY8eOZcOGDRQWFnLqqace\n9MCGzZs3M3/+fDZv3jygx0J8jYk72MSNj65eNwQRCSGEEEfHoQ4+CAaD1NTU8NVXX/HVV19R+eUX\nvPeH3/OAuS3iuZa7Vb6XoadG1fhbW4ilVgO/dgf5UU7k+WJ/6stgrbNpUN/b0fS1J/v96KOPWLFi\nBS+99BJdXV1861vf4sEHH+Scc84JbzN37lxsNhuXXnrpEQdqMBhYvnw5F1xwAcFgkOuvv57x48ez\nYsUKAObNm8esWbNYu3YtZWVlpKSk8Nxzzx30sfFs34kb+35N/FrJ5ap5Nw91aMeEfX8RitiQnMee\n5Dz2JOfRz0L0yI/eZ+X4SaT0dBGoqcLY2kSBLoTdCA6DwvEGqO0IgjmydPHk2qm882eUlpZiX3Yf\nVG7i9GSFX7mD3JSjD3+H3t+qcP5tN8XwnQ6uAbXE6XQ6CgsLmTt3Lj/84Q8pKCiIut2OHTu46aab\n2Lhx41EP9GiJp5Y42Hs499ubN6Nv95CfZKbs1Klxf8h0PJAP2tiTnMee5Dz24jnnhzt1R2dnJ1VV\nVeGpRZw7K/BUfEbl1k08lBV5YOSTbpVzU3TUqlAb0OhIzkBX4CBl1Bjyxx/PJ6/9jbtaP4t43L49\nVPsWiO93h9jcrfFRD6SnpfLNBDs6dUBF3F//+le+/e1vo9dHNkvGm3gr4vo89+tf8YeFt/KAde/f\nYGVqCd+VeeKEEELEQLTWs2dTRnDGTXeSX1xCZWUlVV/tpLXiM7qrd4GrjnRfJ4UGBbux95ysabre\n87H+qS3IXVGmCnnQUsyVv3iE0tJSRo0aRXJy8qFjiDJn6pHOZzccfe0iLpHEaxEn4+OEEELEmqqq\nOJ1OKisreebO21jU9kXENve6VC5O12E3KGQboFmFOlWjNgAuTU8gKxdDQRGpo8ZgHTeBkaNKWb/8\nUW6qfS/iuQbynZZIBdpAfO0xcWLoyblTYy+euzzileQ89iTnsTeYOT/c7k6v18vu3bvDS03lLlor\nPqfbWYnWWIepvRWbXsNmgGRvCHIjy4aQwci7xRP5tKSUvHHHUVJaxnElJXxz5EhsNhs6nS7iMbkp\nyayM1qI2gDHfR3JqsUTdz6WIixMhgynq7bUVnxMMBhOiq1sIIcSRi9bV+MzCCr788guy7UXs3r2b\nul076dhVga+mGp3bRbrPS4FRocAAxQaF4/XQqEKDqlEfgHoNdlkyqLbacblcQGvE644565zD7hHq\nK8Je3KdF7bIEb1EbDNKdGiei/XP+1KVyapKOMlsep//u/yibctrQBSiEECKmuru7cTqd4Va0Vx65\nn6X+3RHb/W+jyqUZegoMYFZ6x6TVqxr1AY0mTY8/Mxe9tZCk4lJyy8YyoqSEESNGMGLECIqLi8Pn\nTB/oeDRxdEl3agKI9qtlxjWn8+iTT3FJTSO2707jjWmzKbnkSv7x3G8S4qS/QgiRCI7kZOyhUIjm\n5ua9XZ3V1TR8VUHbri/x1TrB7SK5uwObQcFqAKtBwdYehJzIr3W/ycLbE04jdeRorGVjKC4p4dQ9\nRVp+fn7U7s5opPVs+JGWuDhxoP78hoYG7rrrLt7/04tcka5jRwB+nLP3H1KOYD1yiTqGYjiTnMee\n5HxwRWu9WqrlcfV9yxgxZhw1NTXU1NRQX1VJR2UFvtrdhJoaMHW0kquEsPUVaAYIAg0quFQNlwpN\nIYVAeja6/EKSi0qo/ORDlgZqImKQA+Diez+XlrgEZrPZeOGFF/jPvHnc+e1ZLM/p7nf/nM4qXlyx\nXIo4IYSIAVVVaWhoCBdnL//8Z/zUW9Vvm4t6XDx/wzWcm6rDalCYaICTld6xaK6gRkOg97LSnMyu\nrHws9hGkjRxNQWkZxcXFTNnTimaz2fqNh964ds0RHywg4pO0xCWQ2y+czvd3vhNx+wM9qSz68985\neeqZQxCVEEIkBr/fT319fbhAq3E6adq1k47qXfjqawi5XZi97eTqNPL0kG9QeL0zyE1RujjvaQqR\nP6IEo81OSvEockeW4SgqwuFw4HA4sNvtpKSkHHaMx9r0G8cCaYk7RmgmS/TbOzy0XzmdJ5JyyPzW\nlVxw653YDnDWDSGESEQHG5emqiqNjY3U1dWFlwbnbtqqvqKnzonaWI/O4ybF5yVPr5BvgDyDwhl6\n8GnQGIQmVaNRhaaQRq0ljbqcfEwFDuorKgBXRDxjpp83KF2cRzL9hohf0hIXJwbSnx9t7MWvTIX0\njJvMu2+WMzXUwTfTdFgU+CyvhILLf8BF188jOzv7iAbeJrp4HkMRryTnsZfIOe87OODVP/2Rfz+x\njNu0vSc9v7/diD8tg2x/D+au9t6WM31vcZZvgFQduFVoDO4pzlSN5pCCLzUDfa4Ni72I9JJSCopH\nhlvPHA4HhYWFmEx7p4SK9rl8L/kseHLFMf8ZG0vxvJ9LS9wxItqRQ9/f05Te3d3NK6+8wgsvv0zV\nxn9ygXetHtwzAAAVEElEQVQXJ/7qHlY/cQ/lyVYCvh5+bGwPP9fKxbv6PacQQgwXmqbR2trar+Ws\nrqaG1upddNZU43PVobU0oe9oI1sXosKn8bP8/l93i9MD3N/UQIlZRxNQZ0qlNisXg7WQZPsIskpK\nKbTbKSws5ITCQgoLC8nLyzvsOTmjziww9Sz5bBVHhbTEHYNaWlp45ZVX+OPLL+Hb+hZZSoh78iPr\n+efGnMVTa98YggiFEIlioK38mqbR0dFBQ0MDDQ0NuFwuGhoaaKxx4q2ppru+lmCzC6WtBVNXB1lK\niFw95BoUcvWQqYf2UG/rWXNQo1kFd1Cjy5xMg0/lvqxgxGs+Yz2ee/64GqvVitFojEU6hDhscu7U\nfUgR15/H4+HWmWcyv60i4r6HmlTK7AWknXwGE2ZfytQZM8nNzR2CKIUQ8Wjj2jX88e7b+KF37wS0\nvwxmkXTGOaSYzXTVOQk01oPHjb7DQ4amkqtXyDXQe6kHowLuYG9B1rynQHMHwWtMgqwcjHk2kgqL\nyCweSYGjiMI9rWaFhYXYbDYsFouce1rEtbjoTm1paeHyyy+nurqakpIS/vSnP5GZmRmx3fr165k/\nfz7BYJA5c+awcOFCAO69915WrlxJXl4eAA888AAXXnhhTN/DYBqs/vzMzExy7EUQpYirCyrYGhsY\n9a9XyP7P//HPAFSb0wmNPg7btHM58cJvcsIJJ2Aw9N+NEmV8XTyPoYhXkvPYO5Kcd3V10dTUhMvl\n2ttiVl9Pq7MKb10NvqYGgi3N1LlcPGjt3/14h76V5Wv/yCWZepqD0Kz2FmXNmka7zkRLWiYVWXmY\n8m0k20eQ7RiB1WbDZrMx3mbDarVitVqxWKIfyBXN7BtuYeXiXcNm6g3Zz2MvUXM+bIq4ZcuWcd55\n53HXXXfx4IMPsmzZMpYtW9Zvm2AwyM0338wbb7yB3W7nlFNOYfbs2YwfPx5FUbj99tu5/fbbh+gd\nxK8DfcD95MUHMWVksnHjRl4s30jHh9sYq3g4ofMd7J9tovk3P2d5QIc7txDDhBMpnDaTkN/P9ud+\nxdx9fnnL+DohYudIfkT1FWVNTU00Njb2Xm904alx0tVQi7+pAbWlGdpbMXg7SQ35ydYrZOohW69Q\nqodJut7uzNYgtAY1WoLQrkRvPWjLL2LHbQux2Wwcb7Vi21OcpaamDkZK5EwDImENm+7UcePG8eab\nb2K1WmloaGD69Ol8/vnn/bbZtGkTS5cuZf369QDhIm/RokUsXbqU1NRU7rjjjoO+jnSnRjeQuYUC\ngQAff/wxmzZtYvOmTVRsfpsM124mmRWOtyiMMik86Q6yMC/yt8GKktP51fp/RbTaCSGOno1r1/Dn\nu+czx1sdvu0JJQ/H7MvIshXSWleDt74Gf2MDgZYmlHYPBm8HKUE/WXqFbD1k6RWy9owx8+5XlLUG\nNVqD0K7oCaZmoGRkYckvIMlmJ9MxAlthIdY9RZnNZmP5/Ju4ruKtiDilG1OIgYuL7lSXy4XVagXA\narXickXOq1NbW0tRUVF43eFwsGXLlvD6U089xQsvvMCUKVP45S9/GbU7VkQ3kLmFjEYjJ510Eied\ndBI33XQTAG63m+3bt/PBBx+w4b1ttP1zNb0nh+mv/d23+bkjGV9mLqYRo0gfPQHrpJMYddzxlJWV\nkZ+fj6Iog/HWhBhUgzl8IBAI0NLSgtvtxu1209zc3Hu9uYn2hnq8rnp87kYCrW60dg/OBhfL8vr/\nH92mNfH4b5/g4hw9fq23KGvZU4z1XXehp86SjpKRjT47jyRrASkFdnKtNvLz88nLy2NCXl74elpa\n2oD+X//7lgWsXLx72HRjCpFoYlrEnXfeeTQ0NETc/otf/KLfuqIoUT8gDvahceONN3LPPfcA8LOf\n/Yw77riD3/72t18z4uFjuPbn5+TkMHPmTGbOnAnAgtkXQpQBxPU6E4HuAA7VRa6nkaLPtmBf9xw9\nGqwPaNRrBnoyclBsDpJGjiZz7HEUjp9IcUkJxcXFpKenh58rVmPuhmvOE1m85TzaHGD7Dx/QNI3O\nzk5aW1vxeDy0trbuXdxuvE0ueppc+NxNBFqbCba1onW2o+/yYgn0kK5XSNdBhh7SdQqj9HCCDrpD\n0BaCtqBGewjagtCihYj2sd6anstfvnkFObaCcCE2bk9R9vnnnzNr1qxB+REl3ZjRxdt+nggSNecx\nLeJef/31A97X141qs9mor68nPz8/Yhu73Y7T6QyvO51OHA4HQL/t58yZw8UXX3zA17r22mspKSkB\negf2T548OfzHLS8vBxh2632GSzwHWh9xxtks3fEJS2gEYFt3iNeSbCx5aQWnnHU2L730ElVVVbQb\njayuqOCL99/F31JHcagHR3cDgfoG8j7exsmpOkwK/LkjRFMQ0i1mejJz+dgPoRY3j2YHws//8A0f\n4Fn2KN+64nu89dZbwyofsn5469u3bz/k9h9s3sTud95E5/dR3e7lzNmXsGDR3YMWn6ZpTJ06lfb2\ndl5//XW8Xi+lpaW0trby1N13clWXE5J0QO/+OLl7Fw9cfTkvZqXT0OrBEPAx1gQZOoVGVSNZB6cn\n6yjRQY8vhBKCE8w62oIa27o1ukOQY4D2IOzwQ8hsJj83B0NGNs3oMWdkMvaEyeTk59Pc3Ex6ejpn\nn302ubm5bL32KrbVfsyUfeIBGD3lZB554smo76+ioiJcwA1G/pTklHDXabx9ng3Wep/hEo+sD6/1\nvutVVVUcyrAZE3fXXXeRk5PDwoULWbZsGR6PJ+LABlVVGTt2LBs2bKCwsJBTTz2VVatWMX78eOrr\n6ynYcyqpxx57jHfffZeXX3454nVkTNzgO5Jz97W0tLBz504qKiqoqqqiuroaV+UufM5d6BvryNcC\nOIwKn/tC3B1lzN3PG1UmJulRk1LRUtNRMrIw5uRhybORVmAn3TGCXMcI8vLzyc3NJSsri4yMDHQ6\n3SHfSyIcaZsoorZ8pZbw3fsfi/i7aJqG1+ulvb09+tLWhrfFTXeLG5+nhUBbK2pHG8HODrSuTuju\nQufrxuD3YSFIik4hRdc7k3+KTiFZgb93BLk5ynkxH2pSmWjR0R7U8ISgR28kmJQKqWno0jMxZOZg\nzM4lJS+fzOwcsrKyyMrKIjs7m5ycHHJzc8nJyRnQPnqo/DybWsJlUfIjhIgPcTFPXEtLC5dddhm7\nd+/uN8VIXV0dc+fOZc2aNQCsW7cuPMXI9ddfz9139/4Cv/rqq9m+fTuKojBy5EhWrFgRHmO3Lyni\n4o+mabjdbqqrq3ly3rXc5vkiYpu7m0KkoZG5p9spc8+Rc5m6vddNCniCe5aQRlsQug0mfOZkgsmp\naKkZ6DKzMWbnYszOpaHVg+c/G/ix3hN+neUGG9MXLOaci75JamoqycnJ/U6xczQNlwJysOLQNA2/\n309XVxddXV14vd7wsu/6vtc3v/g77lNrIp5rUZuZ4rwcdP4eDL4eDKoPkxogSekruvoWhRSl93qy\nDvxa7+D9zhB4Qxpd+617Q9AVgm5Fj2ZJQktKQUlKwZiRgTkzm9odn/ALfXNEPMvtJ7PwN8+TlZVF\nZmbmoO0j0cgJ0IVILHFRxMVKvBZx5Qnan3+4DjZp5/1/fIX6+vr+0yTsc72loZ7uJhdqqxs6PJh6\nuqIWe5l6hUwdPOlWud8WOYv7/U0qZyTr8Gu9RUAAhZDeQMhgRDOaUExmMJnRmcxgtqCYzL2L2YzO\nnITOYkFvSUJvSUaflIwhORm9JRmMJnRGIzq9nqrPd1D32l+4XdcSft3HtBzKLruaCSdNQa/Xo9fr\n0el0Ua9D73kjg8HgIS8Pdt/OTz6iYe3f+LFhbyH7oC+NpClnYrVaCXZ3E+zpJuTrJtTTg+bvQfP5\nwO+DgB8CfpSAH0UNoAuq6INq72UoiCEUxKxTMCu9BbZZgUq/xvEWBYuiYFLApKPf/StbgtxwgJav\nSRYdXg26QtqeIgwCBhOaJQklOQUlJQ19ajrG9ExMmVkkZWaTlplJenr6IRez2Rx1f0yEli/5bIk9\nyXnsxXPO4+LoVCEG4mCTdlosFkaOHMnIkSMH9FyqqtLW1kZLSwutra39Lqvb2+n89eNAa8TjvHoT\n/9an9BYsfh9GLYRJ8WNW/JgUL2Zln8JDp4QLkL7betf33L5fkaLQWxhuaQlye27/f88FiptfPv0w\nx6fvKdIADQgpoO7ZRtln6Vtnz/Pq6f2H73cfoCiRj+lbtrYG+fF+cSw0d/BU+Wucl6nHHwKf1rv4\nNW3P5Z7bQn3XtXDB6w1v23sZ1OnAaEQxW9CZk2g2a9Tk5mBITsGQlIIpJRVTahrmtDTMqWm4XnsF\neqrYX9Kk0/jGk0+TkpJCcnIyKSkppKamDvqUNjJwXwgxlKQlTsSdWHUXDeRUPZqmEQgEwl2CfV1/\n+677fD78fj8+n++Ai6qqvS1ggQCKGmD366+xxNAS8dpLetIonDyFYCjU21oWChEMhghpIdRgKHy7\npmno9HoUnQ6dXo9O0fVf1+1d1+9zufc2A4pOR/W/1nOPLrK78AFzEWfOvQWz2XzQxWQyHfS+wxnv\nBYnR8iWEEIdDulP3IUWcGKihLBiGy7keh0sc+5IxX0KIY8nB6pbD+xkshsz+h6aLwackp/Dd+x/j\nxeNm8lLZmbx43MyYtfjMvuEWVqaW9Lvt2dQSLo7xJKmxjmMg+/mMWRfx6Op1/HL9Rh5dvU4KuK9J\nPltiT3Iee4macxkTJ8RBDORMFoP1ujD0Y62GSxxCCCEiSXeqEEIIIcQwJd2pQgghhBAJRoq4OJGo\n/fnDmeQ89iTnsSc5jz3Jeewlas6liIsTfeeUFLEjOY89yXnsSc5jT3Iee4macyni4oTH4zn0RuKo\nkpzHnuQ89iTnsSc5j71EzbkUcUIIIYQQcUiKuDhRVVU11CEccyTnsSc5jz3JeexJzmMvUXN+zE0x\nMn36dN58882hDkMIIYQQ4pDOPvvsAx6YccwVcUIIIYQQiUC6U4UQQggh4pAUcUIIIYQQcUiKuGFm\n/fr1jBs3jtGjR/Pggw9G3ebWW29l9OjRTJo0iQ8++CDGESaeQ+X8pZdeYtKkSZxwwgmceeaZfPTR\nR0MQZWIZyH4O8O6772IwGPjb3/4Ww+gS00ByXl5ezoknnsjEiROZPn16bANMQIfKeXNzMxdeeCGT\nJ09m4sSJPP/887EPMoH84Ac/wGq1cvzxxx9wm4T7/tTEsKGqqlZaWqpVVlZqfr9fmzRpkrZjx45+\n26xZs0b7xje+oWmapm3evFk77bTThiLUhDGQnL/zzjuax+PRNE3T1q1bJzn/mgaS877tZsyYoV10\n0UXaX/7ylyGINHEMJOetra3ahAkTNKfTqWmapjU1NQ1FqAljIDlfsmSJtmjRIk3TevOdnZ2tBQKB\noQg3Ibz11lva+++/r02cODHq/Yn4/SktccPI1q1bKSsro6SkBKPRyBVXXMHq1av7bfPqq69yzTXX\nAHDaaafh8XhwuVxDEW5CGEjOp06dSkZGBtCb85qamqEINWEMJOcATz31FJdeeil5eXlDEGViGUjO\nX375Zb7zne/gcDgAyM3NHYpQE8ZAcl5QUEB7ezsA7e3t5OTkYDAYhiLchDBt2jSysrIOeH8ifn9K\nETeM1NbWUlRUFF53OBzU1tYechspKo7cQHK+r9/+9rfMmjUrFqElrIHu56tXr+bGG28EQFGUmMaY\naAaS84qKClpaWpgxYwZTpkzh//2//xfrMBPKQHI+d+5cPv30UwoLC5k0aRJPPPFErMM8piTi96eU\n/MPIQL+otP1mhZEvuCN3OLnbuHEjv/vd7/jPf/4ziBElvoHkfP78+SxbtgxFUdA0LWKfF4dnIDkP\nBAK8//77bNiwga6uLqZOncrpp5/O6NGjYxBh4hlIzu+//34mT55MeXk5X331Feeddx4ffvghaWlp\nMYjw2JRo359SxA0jdrsdp9MZXnc6neGujQNtU1NTg91uj1mMiWYgOQf46KOPmDt3LuvXrz9oc704\ntIHk/L333uOKK64Aegd/r1u3DqPRyOzZs2Maa6IYSM6LiorIzc0lKSmJpKQkzjrrLD788EMp4o7Q\nQHL+zjvv8JOf/ASA0tJSRo4cyRdffMGUKVNiGuuxIiG/P4d2SJ7YVyAQ0EaNGqVVVlZqPp/vkAc2\nbNq0KSEGZg6lgeS8urpaKy0t1TZt2jREUSaWgeR8X9dee63217/+NYYRJp6B5Pyzzz7TZs6cqamq\nqnm9Xm3ixInap59+OkQRx7+B5HzBggXavffeq2mapjU0NGh2u11zu91DEW7CqKysHNCBDYny/Skt\nccOIwWBg+fLlXHDBBQSDQa6//nrGjx/PihUrAJg3bx6zZs1i7dq1lJWVkZKSwnPPPTfEUce3geT8\nvvvuo7W1NTw+y2g0snXr1qEMO64NJOfi6BpIzseNG8eFF17ICSecgE6nY+7cuUyYMGGII49fA8n5\n4sWLue6665g0aRKhUIiHHnqI7OzsIY48fn3ve9/jzTffpLm5maKiIpYuXUogEAAS9/tTTrslhBBC\nCBGH5OhUIYQQQog4JEWcEEIIIUQckiJOCCGEECIOSREnhBBCCBGHpIgTQgghhIhDUsQJIYQQQsQh\nKeKEEEIIIeKQFHFCCCGEEHFIijghhBBCiDgkRZwQQhwBr9fLuHHjOO2001BVNXz7P//5T3Q6HU8/\n/fQQRieEOBbIabeEEOIIbd++ndNPP50FCxbwwAMP4HK5mDRpElOnTuWVV14Z6vCEEAlOijghhPga\nHn/8cX784x/zj3/8g4cffphPP/2UDz/8UE5kLoQYdFLECSHE13TRRRexYcMGVFXl9ddfZ8aMGUMd\nkhDiGCBj4oQQ4mv6/ve/j9/vZ9KkSVLACSFiRoo4IYT4GhoaGrjttts4+eST2b59O08++eRQhySE\nOEZIESeEEEdI0zSuueYakpKSeOONN5g/fz4LFy7k448/HurQhBDHABkTJ4QQR+iRRx5h0aJFbNy4\nkWnTphEIBDj99NPx+Xxs27YNi8Uy1CEKIRKYtMQJIcQReP/99/npT3/K4sWLmTZtGgBGo5FVq1ZR\nVVXFHXfcMcQRCiESnbTECSGEEELEIWmJE0IIIYSIQ1LECSGEEELEISnihBBCCCHikBRxQgghhBBx\nSIo4IYQQQog4JEWcEEIIIUQckiJOCCGEECIOSREnhBBCCBGHpIgTQgghhIhD/x8YNYW4nOwjTwAA\nAABJRU5ErkJggg==\n", "text": [ "" ] } ], "prompt_number": 5 }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ "Free-stream conditions" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The airfoil is immersed in a free-stream ($U_\\infty$,$\\alpha$) where $U_\\infty$ and $\\alpha$ are the velocity magnitude and angle of attack, respectively. Like before, we create a class for the free stream, even though we will only have one object that uses this class. It makes it easier to pass the free stream to other functions later on." ] }, { "cell_type": "code", "collapsed": false, "input": [ "class Freestream:\n", " \"\"\"Freestream conditions.\"\"\"\n", " def __init__(self, u_inf=1.0, alpha=0.0):\n", " \"\"\"Sets the freestream conditions.\n", " \n", " Arguments\n", " ---------\n", " u_inf -- Farfield speed (default 1.0).\n", " alpha -- Angle of attack in degrees (default 0.0).\n", " \"\"\"\n", " self.u_inf = u_inf\n", " self.alpha = alpha*pi/180 # degrees --> radians" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 6 }, { "cell_type": "code", "collapsed": false, "input": [ "# defines and creates the object freestream\n", "u_inf = 1.0 # freestream spee\n", "alpha = 0.0 # angle of attack (in degrees)\n", "freestream = Freestream(u_inf, alpha) # instantiation of the object freestream" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 7 }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ "Flow tangency boundary condition" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A constant vortex strength $\\gamma$ will be added to each panel (all panels have the same, constant vortex-sheet strength). Thus, using the principle of superposition, the velocity potential becomes:\n", "\n", "$$\n", "\\begin{align*}\n", "\\phi\\left(x_{c_i},y_{c_i}\\right) &= V_\\infty x_{c_i} \\cos \\alpha + V_\\infty y_{c_i} \\sin \\alpha \\\\\n", "&+ \\sum_{j=1}^N \\frac{\\sigma_j}{2\\pi} \\int_j \\ln \\left(\\sqrt{(x_{c_i}-x_j(s_j))^2+(y_{c_i}-y_j(s_j))^2} \\right) {\\rm d}s_j \\\\\n", "&- \\sum_{j=1}^N \\frac{\\gamma}{2\\pi} \\int_j \\tan^{-1} \\left(\\frac{y_{c_i}-y_j(s_j)}{x_{c_i}-x_j(s_j)}\\right) {\\rm d}s_j\n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The flow tangency boundary condition is applied at every panel center:\n", "\n", "$$0 = \\underline{V}\\cdot\\underline{n}_i = \\frac{\\partial}{\\partial n_i} \\left\\{ \\phi\\left(x_{c_i},y_{c_i}\\right) \\right\\}$$\n", "\n", "i.e.\n", "\n", "$$\n", "\\begin{align*}\n", "0 &= V_\\infty \\cos \\left(\\alpha-\\beta_i\\right) + \\frac{\\sigma_i}{2} \\\\\n", "&+ \\sum_{j=1,j\\neq i}^N \\frac{\\sigma_j}{2\\pi} \\int_j \\frac{\\partial}{\\partial n_i} \\ln \\left(\\sqrt{(x_{c_i}-x_j(s_j))^2+(y_{c_i}-y_j(s_j))^2} \\right) {\\rm d}s_j \\\\\n", "&- \\sum_{j=1,j\\neq i}^N \\frac{\\gamma}{2\\pi} \\int_j \\frac{\\partial}{\\partial n_i} \\tan^{-1} \\left(\\frac{y_{c_i}-y_j(s_j)}{x_{c_i}-x_j(s_j)}\\right) {\\rm d}s_j\n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We already worked the first integral in the previous lesson:\n", "\n", "$$\\frac{\\partial}{\\partial n_i} \\ln \\left(\\sqrt{(x_{c_i}-x_j(s_j))^2+(y_{c_i}-y_j(s_j))^2} \\right) = \\frac{\\left(x_{c_i}-x_j\\right)\\frac{\\partial x_{c_i}}{\\partial n_i} + \\left(y_{c_i}-y_j\\right)\\frac{\\partial y_{c_i}}{\\partial n_i}}{\\left(x_{c_i}-x_j\\right)^2 + \\left(x_{c_i}-x_j\\right)^2}$$\n", "\n", "where $\\frac{\\partial x_{c_i}}{\\partial n_i} = \\cos \\beta_i$ and $\\frac{\\partial y_{c_i}}{\\partial n_i} = \\sin \\beta_i$, and:\n", "\n", "$$x_j(s_j) = x_{b_j} - s_j \\sin \\beta_j$$\n", "\n", "$$y_j(s_j) = y_{b_j} + s_j \\cos \\beta_j$$\n", "\n", "We now need to derive the last integral of the boundary equation:\n", "\n", "$$\\frac{\\partial}{\\partial n_i} \\tan^{-1} \\left(\\frac{y_{c_i}-y_j(s_j)}{x_{c_i}-x_j(s_j)}\\right)= \\frac{\\left(x_{c_i}-x_j\\right)\\frac{\\partial y_{c_i}}{\\partial n_i} - \\left(y_{c_i}-y_j\\right)\\frac{\\partial x_{c_i}}{\\partial n_i}}{\\left(x_{c_i}-x_j\\right)^2 + \\left(y_{c_i}-y_j\\right)^2}$$\n", "\n", "where $\\frac{\\partial x_{c_i}}{\\partial n_i} = \\cos \\beta_i$ and $\\frac{\\partial y_{c_i}}{\\partial n_i} = \\sin \\beta_i$." ] }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ "Enforcing the Kutta-condition" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To enforce the *Kutta-condition*, we state that the pressure coefficient on the fisrt panel must be equal to that on the last panel:\n", "\n", "$$C_{p_1} = C_{p_{N}}$$\n", "\n", "Using the definition of the pressure coefficient $C_p = 1-\\left(\\frac{V}{U_\\infty}\\right)^2$, the Kutta-condition implies that the magnitude of the velocity at the first panel center must equal the magnitude of the last panel center:\n", "\n", "$$V_1^2 = V_N^2$$\n", "\n", "Since the flow tangency condition requires that $V_{n_1} = V_{n_N} = 0$, we end up with the following *Kutta-condition*:\n", "\n", "$$V_{t_1} = - V_{t_N}$$\n", "\n", "(the minus sign comes from the reference axis we chose for the normal and tangential vectors).\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Therefore, we need to evaluate the tangential velocity at the first and last panels. Let's derive it for every panel, since it will be useful to compute the pressure coefficient.\n", "\n", "$$V_{t_i} = \\frac{\\partial}{\\partial t_i} \\left(\\phi\\left(x_{c_i},y_{c_i}\\right)\\right)$$\n", "\n", "i.e.,\n", "\n", "$$\n", "\\begin{align*}\n", "V_{t_i} &= V_\\infty \\sin \\left(\\alpha-\\beta_i\\right) \\\\\n", "&+ \\sum_{j=1,j\\neq i}^N \\frac{\\sigma_j}{2\\pi} \\int_j \\frac{\\partial}{\\partial t_i} \\ln \\left(\\sqrt{(x_{c_i}-x_j(s_j))^2+(y_{c_i}-y_j(s_j))^2} \\right) {\\rm d}s_j \\\\\n", "&- \\sum_{j=1,j\\neq i}^N \\frac{\\gamma}{2\\pi} \\int_j \\frac{\\partial}{\\partial t_i} \\tan^{-1} \\left(\\frac{y_{c_i}-y_j(s_j)}{x_{c_i}-x_j(s_j)}\\right) {\\rm d}s_j\n", "\\end{align*}\n", "$$\n", "\n", "\n", "which gives\n", "\n", "$$\n", "\\begin{align*}\n", "V_{t_i} &= V_\\infty \\sin \\left(\\alpha-\\beta_i\\right) \\\\\n", "&+ \\sum_{j=1,j\\neq i}^N \\frac{\\sigma_j}{2\\pi} \\int_j \\frac{\\left(x_{c_i}-x_j\\right)\\frac{\\partial x_{c_i}}{\\partial t_i} + \\left(y_{c_i}-y_j\\right)\\frac{\\partial y_{c_i}}{\\partial t_i}}{\\left(x_{c_i}-x_j\\right)^2 + \\left(x_{c_i}-x_j\\right)^2} {\\rm d}s_j \\\\\n", "&- \\sum_{j=1,j\\neq i}^N \\frac{\\gamma}{2\\pi} \\int_j \\frac{\\left(x_{c_i}-x_j\\right)\\frac{\\partial y_{c_i}}{\\partial t_i} - \\left(y_{c_i}-y_j\\right)\\frac{\\partial x_{c_i}}{\\partial t_i}}{\\left(x_{c_i}-x_j\\right)^2 + \\left(x_{c_i}-x_j\\right)^2} {\\rm d}s_j \n", "\\end{align*}\n", "$$\n", "\n", "where $\\frac{\\partial x_{c_i}}{\\partial t_i} = -\\sin \\beta_i$ and $\\frac{\\partial y_{c_i}}{\\partial t_i} = \\cos \\beta_i$" ] }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ "Building the linear system" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here, we build and solve the linear system of equations of the form\n", "\n", "$$[A][\\sigma,\\gamma] = [b]$$\n", "\n", "where the $N+1 \\times N+1$ matrix $[A]$ contains three blocks: an $N \\times N$ source matrix (the same one of Lesson 10), an $N \\times 1$ vortex array to store the weight of the variable $\\gamma$ at each panel, and a $1 \\times N+1$ Kutta array that repesents our Kutta-condition." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We are going to re-use the function `integral()` from Lesson 10 to compute the different integrals with the NumPy function `integrate.quad()`:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "def integral(x, y, panel, dxdz, dydz):\n", " \"\"\"Evaluates the contribution of a panel at one point.\n", " \n", " Arguments\n", " ---------\n", " x, y -- Cartesian coordinates of the point.\n", " panel -- panel which contribution is evaluated.\n", " dxdz -- derivative of x in the z-direction.\n", " dydz -- derivative of y in the z-direction.\n", " \n", " Returns\n", " -------\n", " Integral over the panel of the influence at one point.\n", " \"\"\"\n", " def func(s):\n", " return ( ((x - (panel.xa - sin(panel.beta)*s))*dxdz \n", " + (y - (panel.ya + cos(panel.beta)*s))*dydz)\n", " / ((x - (panel.xa - sin(panel.beta)*s))**2 \n", " + (y - (panel.ya + cos(panel.beta)*s))**2) )\n", " return integrate.quad(lambda s:func(s), 0., panel.length)[0]" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 8 }, { "cell_type": "markdown", "metadata": {}, "source": [ "We first build the source matrix:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "def source_matrix(panels):\n", " \"\"\"Builds the source matrix.\n", " \n", " Arguments\n", " ---------\n", " panels -- array of panels.\n", " \n", " Returns\n", " -------\n", " A -- NxN matrix (N is the number of panels).\n", " \"\"\"\n", " N = len(panels)\n", " A = np.empty((N, N), dtype=float)\n", " np.fill_diagonal(A, 0.5)\n", " \n", " for i, p_i in enumerate(panels):\n", " for j, p_j in enumerate(panels):\n", " if i != j:\n", " A[i,j] = 0.5/pi*integral(p_i.xc, p_i.yc, p_j, cos(p_i.beta), sin(p_i.beta))\n", " \n", " return A" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 9 }, { "cell_type": "markdown", "metadata": {}, "source": [ "then the vortex array:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "def vortex_array(panels):\n", " \"\"\"Builds the vortex array.\n", " \n", " Arguments\n", " ---------\n", " panels - array of panels.\n", " \n", " Returns\n", " -------\n", " a -- 1D array (Nx1, N is the number of panels).\n", " \"\"\"\n", " a = np.zeros(len(panels), dtype=float)\n", " \n", " for i, p_i in enumerate(panels):\n", " for j, p_j in enumerate(panels):\n", " if i != j:\n", " a[i] -= 0.5/pi*integral(p_i.xc, p_i.yc, p_j, +sin(p_i.beta), -cos(p_i.beta))\n", " \n", " return a" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 10 }, { "cell_type": "markdown", "metadata": {}, "source": [ "and finally the Kutta array:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "def kutta_array(panels):\n", " \"\"\"Builds the Kutta-condition array.\n", " \n", " Arguments\n", " ---------\n", " panels -- array of panels.\n", " \n", " Returns\n", " -------\n", " a -- 1D array (Nx1, N is the number of panels).\n", " \"\"\"\n", " N = len(panels)\n", " a = np.zeros(N+1, dtype=float)\n", " \n", " a[0] = 0.5/pi*integral(panels[N-1].xc, panels[N-1].yc, panels[0], \n", " -sin(panels[N-1].beta), +cos(panels[N-1].beta))\n", " a[N-1] = 0.5/pi*integral(panels[0].xc, panels[0].yc, panels[N-1], \n", " -sin(panels[0].beta), +cos(panels[0].beta))\n", " \n", " for i, panel in enumerate(panels[1:N-1]):\n", " a[i] = 0.5/pi*(integral(panels[0].xc, panels[0].yc, panel, \n", " -sin(panels[0].beta), +cos(panels[0].beta))\n", " + integral(panels[N-1].xc, panels[N-1].yc, panel, \n", " -sin(panels[N-1].beta), +cos(panels[N-1].beta)) )\n", "\n", " a[N] -= 0.5/pi*(integral(panels[0].xc, panels[0].yc, panel, \n", " +cos(panels[0].beta), +sin(panels[0].beta))\n", " + integral(panels[N-1].xc, panels[N-1].yc, panel, \n", " +cos(panels[N-1].beta), +sin(panels[N-1].beta)) )\n", " \n", " return a" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 11 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now that the three blocks have be defined, we can assemble the matrix $[A]$:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "def build_matrix(panels):\n", " \"\"\"Builds the matrix of the linear system.\n", " \n", " Arguments\n", " ---------\n", " panels -- array of panels.\n", " \n", " Returns\n", " -------\n", " A -- (N+1)x(N+1) matrix (N is the number of panels).\n", " \"\"\"\n", " N = len(panels)\n", " A = np.empty((N+1, N+1), dtype=float)\n", " \n", " AS = source_matrix(panels)\n", " av = vortex_array(panels)\n", " ak = kutta_array(panels)\n", " \n", " A[0:N,0:N], A[0:N,N], A[N,:] = AS[:,:], av[:], ak[:]\n", " \n", " return A" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 12 }, { "cell_type": "markdown", "metadata": {}, "source": [ "On the right hand-side, we store the free-stream conditions that do not depend on the unknowns strengths:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "def build_rhs(panels, freestream):\n", " \"\"\"Builds the RHS of the linear system.\n", " \n", " Arguments\n", " ---------\n", " panels -- array of panels.\n", " freestream -- farfield conditions.\n", " \n", " Returns\n", " -------\n", " b -- 1D array ((N+1)x1, N is the number of panels).\n", " \"\"\"\n", " N = len(panels)\n", " b = np.empty(N+1,dtype=float)\n", " \n", " for i, panel in enumerate(panels):\n", " b[i] = - freestream.u_inf * cos(freestream.alpha - panel.beta)\n", " b[N] = -freestream.u_inf*( sin(freestream.alpha-panels[0].beta)\n", " +sin(freestream.alpha-panels[N-1].beta) )\n", " \n", " return b" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 13 }, { "cell_type": "markdown", "metadata": {}, "source": [ "All in all, the system has been defined with the functions `source_matrix()`, `vortex_array()`, `kutta_array()` and `build_rhs()`, which we have to call to build the complete system:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "A = build_matrix(panels) # calculates the singularity matrix\n", "b = build_rhs(panels, freestream) # calculates the freestream RHS" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 14 }, { "cell_type": "markdown", "metadata": {}, "source": [ "The linear system is then solved using the NumPy function `linalg.solve()`and we store the results in the attribute `sigma` of each object. We also create a variable `gamma` to store the value of the constant vortex strength. " ] }, { "cell_type": "code", "collapsed": false, "input": [ "# solves the linear system\n", "variables = np.linalg.solve(A, b)\n", "\n", "for i, panel in enumerate(panels):\n", "\tpanel.sigma = variables[i]\n", "gamma = variables[-1]" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 15 }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ "Surface pressure coefficient" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The pressure coefficient at the $i$-th panel center is:\n", "\n", "$$C_{p_i} = 1 - \\left(\\frac{V_{t_i}}{U_\\infty}\\right)^2$$\n", "\n", "So, we have to compute the tangential velocity at each panel center using the function `get_tangential_velocity()`:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "def get_tangential_velocity(panels, freestream, gamma):\n", " \"\"\"Computes the tangential velocity on the surface.\n", " \n", " Arguments\n", " ---------\n", " panels -- array of panels.\n", " freestream -- farfield conditions.\n", " gamma -- circulation density.\n", " \"\"\"\n", " N = len(panels)\n", " A = np.empty((N, N+1), dtype=float)\n", " np.fill_diagonal(A, 0.0)\n", " \n", " for i, p_i in enumerate(panels):\n", " for j, p_j in enumerate(panels):\n", " if i != j:\n", " A[i,j] = 0.5/pi*integral(p_i.xc, p_i.yc, p_j, -sin(p_i.beta), +cos(p_i.beta))\n", " A[i,N] -= 0.5/pi*integral(p_i.xc, p_i.yc, p_j, +cos(p_i.beta), +sin(p_i.beta))\n", "\n", " b = freestream.u_inf * np.sin([freestream.alpha - panel.beta for panel in panels])\n", " \n", " var = np.append([panel.sigma for panel in panels], gamma)\n", " \n", " vt = np.dot(A, var) + b\n", " for i, panel in enumerate(panels):\n", " panel.vt = vt[i]" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 16 }, { "cell_type": "code", "collapsed": false, "input": [ "# computes the tangential velocity at each panel center.\n", "get_tangential_velocity(panels, freestream, gamma)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 17 }, { "cell_type": "markdown", "metadata": {}, "source": [ "And we define a function `get_pressure_coefficient()` to compute the surface pressure coefficient:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "def get_pressure_coefficient(panels, freestream):\n", " \"\"\"Computes the surface pressure coefficients.\n", " \n", " Arguments\n", " ---------\n", " panels -- array of panels.\n", " freestream -- farfield conditions.\n", " \"\"\"\n", " for panel in panels:\n", " panel.cp = 1.0 - (panel.vt/freestream.u_inf)**2" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 18 }, { "cell_type": "code", "collapsed": false, "input": [ "# computes surface pressure coefficient\n", "get_pressure_coefficient(panels, freestream)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 19 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Time to plot the result!" ] }, { "cell_type": "code", "collapsed": false, "input": [ "# plots the surface pressure coefficient\n", "val_x, val_y = 0.1, 0.2\n", "x_min, x_max = min( panel.xa for panel in panels ), max( panel.xa for panel in panels )\n", "cp_min, cp_max = min( panel.cp for panel in panels ), max( panel.cp for panel in panels )\n", "x_start, x_end = x_min-val_x*(x_max-x_min), x_max+val_x*(x_max-x_min)\n", "y_start, y_end = cp_min-val_y*(cp_max-cp_min), cp_max+val_y*(cp_max-cp_min)\n", "\n", "plt.figure(figsize=(10, 6))\n", "plt.grid(True)\n", "plt.xlabel('x', fontsize=16)\n", "plt.ylabel('$C_p$', fontsize=16)\n", "plt.plot([panel.xc for panel in panels if panel.loc == 'extrados'], \n", " [panel.cp for panel in panels if panel.loc == 'extrados'], \n", " color='r', linestyle='-', linewidth=2, marker='o', markersize=6)\n", "plt.plot([panel.xc for panel in panels if panel.loc == 'intrados'], \n", " [panel.cp for panel in panels if panel.loc == 'intrados'], \n", " color='b', linestyle='-', linewidth=1, marker='o', markersize=6)\n", "plt.legend(['extrados', 'intrados'], loc='best', prop={'size':14})\n", "plt.xlim(x_start, x_end)\n", "plt.ylim(y_start, y_end)\n", "plt.gca().invert_yaxis()\n", "plt.title('Number of panels : %d' % N);" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAnEAAAGOCAYAAAD4oVVpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XlcVPX6wPHPsIOyq6ggkrjihuKu4KABKZm2YaXX5aaZ\nFWr3mlbupaXXcsP6VZp1SzNbrFuiiBvgmoo7gju4ILggIvt2fn9Mjo6gouIMMzzv1+u8bufMOWee\neZgrD+f7Pc9RKYqiIIQQQgghjIqZoQMQQgghhBAPToo4IYQQQggjJEWcEEIIIYQRkiJOCCGEEMII\nSREnhBBCCGGEpIgTQgghhDBCUsQJIYzKsGHDmDJlisHef/jw4bi4uNClSxeDxXCnmJgYGjRoYOgw\nhBB6JkWcEOKReHl54ebmRm5urnbb0qVLCQwMfCzvp1KpUKlUj+Xc97N161Y2btxIamoqu3btMkgM\nhpCRkUHt2rXx9/fX2X7gwAH8/PyoUaMGHTp04ODBgwaKUIjqSYo4IcQjKy0tZeHChXp7v8rqUV5a\nWvpA+6ekpODl5YWNjU2lvL+xmDhxIj4+PjrFc2FhIf3792fIkCFkZmYydOhQ+vfvT1FRkQEjFaJ6\nkSJOCPFIVCoV48eP55NPPuH69etlXk9OTsbMzEynYFKr1Xz99dcAfPvtt3Tv3p1//etfODs707hx\nY3bs2ME333yDp6cnbm5ufPfddzrnvHLlCsHBwTg4OKBWqzl79qz2taSkJIKCgnB1daV58+b8/PPP\n2teGDRvG6NGj6du3LzVr1iQmJqZMvKmpqTzzzDO4urrSpEkTli5dCsDXX3/NyJEj2blzJ/b29syY\nMaPMsTc/S3h4OE5OTrRo0YLNmzdrX//mm2/w8fHBwcEBb29vvvrqK+1rMTExeHh4MG/ePNzc3Khf\nvz7ffvut9vWCggLGjx9Pw4YNqVu3LqNHjyY/P7/cn8mcOXPw8PDAwcGB5s2b68TwoHbs2EFCQgLD\nhw/XKZ5jYmIoKSlh7NixWFpaEh4ejqIoj/ReQogHI0WcEOKRdejQAbVazSeffFKh/e8cEt29ezdt\n27YlIyODl19+mbCwMPbt28epU6dYvnw5b731lna4VlEUVqxYwdSpU7ly5Qq+vr4MGjQIgJycHIKC\nghg8eDCXL1/mxx9/5I033iAxMVH7XitXrmTKlClkZ2fTvXv3MrG99NJLeHp6cvHiRX755Rfef/99\ntmzZwquvvsoXX3xB165duXHjBtOmTSv3s+3evZvGjRtz9epVZsyYwXPPPce1a9cAcHNzIzIykqys\nLL755hvefvtt9u/frz02PT2drKwsUlNT+frrr3nzzTe1hfG7777LyZMnOXjwICdPnuTChQt88MEH\nZd7/2LFjfPbZZ+zdu5esrCyio6Px8vIqN9YffviBtm3b3vXnVFJSQnh4OJ999lmZ1xISEmjTpo3O\ntrZt25KQkHDX8wkhKpcUcUKIR6ZSqfjggw+IiIjgypUrD3z8E088wdChQ1GpVISFhZGamsrUqVOx\ntLQkKCgIKysrTp48qd3/6aefpkePHlhZWTFr1ix27tzJ+fPnWbNmjfZcZmZm+Pr68txzz+lcjRsw\nYABdu3YFwNraWieOc+fOsWPHDubMmYOVlRVt27ZlxIgR2iuBFRnGrVOnDmPHjsXc3JywsDCaNWtG\nZGQkAH379uWJJ54AICAggODgYLZu3ao91tLSkqlTp2Jubk6fPn2oWbMmx44dQ1EUlixZwrx583By\ncqJmzZq89957/Pjjj2Xe39zcnIKCAhISEigqKsLT05NGjRqVG+srr7xyz3lsixYtokuXLrRr167M\na9nZ2Tg6Oupsc3Bw4MaNG/fNkRCiclgYOgAhhGlo2bIlTz/9NLNnz6ZFixYPdKybm5v2v21tbQGo\nXbu2zrbs7GxAUzB6eHhoX6tRowYuLi6kpqaSkpLCX3/9hbOzs/b14uJihgwZUu6xd0pNTcXFxYUa\nNWpot3l6erJ3794KfxZ3d3ed9YYNG3Lx4kUA1q1bx4wZMzhx4gSlpaXk5ubqXM1ydXXFzOzW39Z2\ndnZkZ2dz+fJlcnNz8fPz076mKEq5c/oaN27MggULmD59OgkJCYSEhDBv3jzq1atX4c8AmlxEREQQ\nHx9f7uv29vZkZWXpbLt+/ToODg4P9D5CiIcnV+KEEJVmxowZLFmyhAsXLmi33SyIbr97NS0t7aHf\nQ1EUzp07p13Pzs4mIyMDd3d3PD096dmzJ9euXdMuN27cKHc4sDz169cnIyNDWzACnD179p6F351u\n/+yguRmifv36FBQU8PzzzzNhwgQuXbrEtWvX6Nu3b4Wu7tWqVQtbW1uOHj2q/VyZmZlliqibXn75\nZbZu3UpKSgoqlYqJEydWOP6bdu/ezcWLF/Hx8aFevXqMGzeO3bt3U79+fRRFwcfHh0OHDukcc+jQ\nIVq2bPnA7yWEeDhSxAkhKo23tzcDBw7UuVO1du3auLu78/3331NSUsKyZcs4derUI73P2rVr2b59\nO4WFhUyZMoWuXbvi7u5OaGgox48fZ/ny5RQVFVFUVMSePXtISkoC7j8c2qBBA7p168Z7771HQUEB\nhw4dYtmyZQwePLjCsV26dIlFixZRVFTEzz//TFJSEn379qWwsJDCwkJq1aqFmZkZ69atIzo6ukLn\nNDMzY+TIkYwbN47Lly8DmmKxvOOPHz/O5s2bKSgowNraGhsbG8zNzSsc/019+/YlJSWFgwcPcvDg\nQT744APatWvHgQMHUKlUBAYGYm5uzqJFiygoKGDRokWYmZnRq1evB34vIcTDkSJOCFGppk6dSm5u\nrs6NC0uWLGHu3LnUqlWLo0eP6txQUF7ft3v1gVOpVAwaNIgZM2bg6urK/v37Wb58OaAZ4ouOjubH\nH3/E3d2devXq8d5771FYWHjX97rTypUrSU5Opn79+jz33HN88MEH2sKkIsd37tyZEydOULt2baZM\nmcKvv/6Ks7Mz9vb2LFq0iLCwMFxcXFi5ciX9+/ev8OeeM2cOjRs3pkuXLjg6OhIUFMTx48fLHFtQ\nUMB7771H7dq1qVevHleuXOHjjz8u95wrVqygVatW5b5mZWVFnTp1tIujo6N2G2jm7/3+++989913\nODs789133/H7779jYSGzdITQF5VSWQ2XhBCimvv222/5+uuvdW5WEEKIx0WuxAkhhBBCGCEp4oQQ\nopIY8pFgQojqR4ZThRBCCCGMkFyJE0IIIYQwQtXuNiK1Wk1sbKyhwxBCCCGEuK+ePXuW+5xnqIZX\n4mJjY1EUxeiWadOmGTyG6rZIziXn1WGRnEvOq8NizDm/14WnalfEGavk5GRDh1DtSM71T3Kuf5Jz\n/ZOc65+p5lyKOCGEEEIIIyRFnJEYNmyYoUOodiTn+ic51z/Juf5JzvXPVHNe5VuMZGRkMHDgQFJS\nUvDy8uKnn37CycmpzH6ZmZmMGDGChIQEVCoVy5Yto0uXLmX2U6lUVPGPLIQQQggB3LtuqfJX4mbP\nnq19RmDv3r2ZPXt2ufuNHTuWvn37kpiYyKFDh2jRooWeI3287nZninh8JOf6JznXP8m5/knO9c9U\nc17lW4z88ccf2jszhg4dilqtLlPIXb9+na1bt/Lf//4XAAsLCxwdHfUeqxBCCGFILi4uXLt2zdBh\niAfk7OxMRkbGAx9X5YdTnZ2dtV9IRVHK/YIeOHCAUaNG4ePjw8GDB/Hz82PhwoXY2dmVOZ8Mpwoh\nhDBV8jvOON3r51blh1ODgoJo3bp1meWPP/7Q2e9uzyUsLi5m3759vPHGG+zbt48aNWrcddhVCCGE\nEMIUVInh1A0bNtz1NTc3N9LS0qhbty4XL16kTp06Zfbx8PDAw8ODjh07AvDCCy/cs4gbNmwYXl5e\nADg5OeHr64tarQZujZtXtfWb26pKPNVh/c7cGzqe6rC+YMECo/j/oymtHzhwgHHjxlWZeKrD+s1t\nj/P8wvjc/vOLiYmpUG+7Kj+cOmHCBFxdXZk4cSKzZ88mMzOz3AItICCApUuX0rRpU6ZPn05eXh5z\n5swps5+xXmqOiYnR/h9V6IfkXP8k5/onOde/x5lzY/0dV9097HBqlS/iMjIyCAsL4+zZszotRlJT\nUxk5ciSRkZEAHDx4kBEjRlBYWIi3tzfffPNNuTc3yBdcCCGEqZLfccbJZIu4yiZfcCGEEKZKfsdV\nvr1799KpUyeSk5Px9PR8LO9h1Dc2iPuTuQ76JznXP8m5/knO9U9yXnmSk5MxMzNj3759hg7FIKSI\nE0IIIYRRu9/Vx8LCQj1Fol9SxBkJmXisf5Jz/ZOc65/kXP8MlfO4yEgmh4QwXa1mckgIcX/PKa8K\n5/vPf/5D48aNsbOzo02bNqxYsQKAV199lVatWpGfnw9ASUkJ/v7+PPPMMwA0atQIgI4dO2JmZkav\nXr0ATReKfv36MWfOHDw8PLTDoMuXL6djx444ODjg5uZGWFgYqampOrFERUXRvHlzbG1tCQgI4Pjx\n42XiXb16Na1bt8bGxgZPT08++uijMq+3adMGOzs7XF1dUavVXLp06aHzc1dKNVMNP7IQQohq4m6/\n42LXrFHe9/ZWFNAu73t7K7Fr1jzU+1Tm+d5//32lefPmyvr165Xk5GTlhx9+UGrUqKFERkYqOTk5\nStOmTZU333xTURRFmTFjhlKvXj3l8uXLiqIoyp49exSVSqVER0cr6enpyrVr1xRFUZShQ4cq9vb2\nyuDBg5WEhATlyJEjiqIoyrJly5R169YpZ86cUXbv3q0EBgYqAQEB2ljOnj2rWFtbK2PGjFGOHTum\n/PTTT4q7u7tiZmampKSkKIqiKHv37lXMzc2V6dOnKydOnFBWrFih1KxZU4mIiFAURVEuXryoWFpa\nKvPmzVNSUlKUI0eOKF9//bWSnp5+1xzcqza552sVSbApMdYibsuWLYYOodqRnOuf5Fz/JOf69zhz\nfrffcZOCg3UKrpvL5HK2VWSZdJftk0NCHije7OxsxdbWVtm2bZvO9rFjxyp9+/ZVFEVTqFlZWSlT\npkxRLC0tlaioKO1+Z86cUVQqlRIfH69z/NChQ5U6deoohYWF93z/xMRERaVSKRcuXFAURVHee+89\npVmzZjr7zJw5U1GpVNoi7pVXXlF69+6ts8/06dMVDw8PRVEUJT4+Xmf/injYIk6GU4UQQggTZ1FQ\nUO5284c93122m/897FlRR48eJT8/n5CQEOzt7bXLF198wenTpwHo0KEDkyZNYubMmYwaNYqQkJAK\nnbtVq1ZYWlrqbNu3bx/9+/fHy8sLBwcH7UMCzp49C0BiYiJdunTROebO9aSkJLp3766zrXv37ly4\ncIHs7Gx8fX158sknadWqFS+88AJffPEFV65cqXhSHoAUcUZC5q3on+Rc/yTn+ic51z9D5LzY2rrc\n7SUhIQ91La44OLj889nYPFBcpaWlAKxZs4aDBw9ql6NHjxIdHQ1oblrYunUr5ubmnDx5ssLnvvP5\n6Tk5OYSEhFCzZk2WL1/O3r17iYqKAm7d+FDRFi1320elUmFmZkZ0dDTR0dG0adOGr7/+miZNmnDo\n0KEKx15RUsQJIYQQJi54zBgmeXvrbHvf25ug8HCDns/Hxwdra2uSk5Np1KiRztKgQQMA5s2bx4ED\nB9i6dSu7du0iIiJCe7yVlRWgueHhfpKSkrh69SofffQRPXr0oGnTpqSnp+vs06JFC/766y+dbbt2\n7Sqzz/bt23W2bdu2jQYNGlCjRg3tti5dujB16lT27NlD/fr1WbVqVQUy8mCqxLNTxf3Jo3H0T3Ku\nf5Jz/ZOc658hch4QGgrAlIgIzPPzKbGx4anwcO12Q53P3t6e8ePHM378eBRFwd/fn+zsbHbt2oW5\nuTmdOnVi8uTJrFy5ki5duvD555/z6quv0rt3b3x8fKhTpw62trZERUXh6emJra0tDg4O5b6Xp6cn\n1tbWRERE8MYbb5CYmMiUKVN09nn99df59NNPGTduHKNHj+bw4cN8+eWXOvv8+9//pmPHjsyYMYOX\nX36ZPXv2MG/ePD7++GNAU/Rt3LiRp556ijp16rB//37OnTtHy5YtHyg3FVLhWXcmwlg/skw+1j/J\nuf5JzvVPcq5/hrixoaqLiIhQfHx8FGtra6V27dpKcHCw8ueffyqtWrVSXn31VZ19//GPfyi+vr7a\nmxaWLl2qeHp6Kubm5kpgYKCiKIoybNgwpV+/fmXeZ9WqVYq3t7diY2OjdO7cWVm/fr1iZmamxMbG\naveJjIxUmjVrptjY2Cg9evRQVqxYoXN3qqIoyurVq5XWrVsrVlZWiqenp/LRRx9pX0tMTFT69Omj\nuLm5KdbW1kqTJk2UuXPn3vPz3+vndq/X5LFbQgghhImQ33HGSR67JYQQQghRjUgRZyTkWXv6JznX\nP8m5/knO9U9yLiqLFHFCCCGEEEZI5sQJIYQQJkJ+xxknmRMnhBBCCFGNSBFnJGQOhf5JzvVPcq5/\nknP9k5yLyiJFnBBCCCGEEZI5cUIIIYSJkN9xxknmxAkhhBBCVCNSxBkJmUOhf5Jz/ZOc65/kXP8k\n57cMGzaMfv36GToMHTVr1uS///2vocOoECnihBBCCGEQERERrFixosL7e3l58emnnz7GiDTDlyqV\n6rG+R2WxMHQAomLUarWhQ6h2JOf6JznXP8m5/knOb7G3t3+g/StSXJWWlgJgZmb616lM/xMKIYQQ\ngsjIOEJCJqNWTyckZDKRkXEGP9/tw6lqtZo333yT999/n9q1a+Pm5sY777yjndSvVqtJSUnhnXfe\nwczMDHNzcwC+/fZb7O3tWbduHa1atcLa2pqkpCT27NlDcHAwtWvXxtHREX9/f3bt2qXz/idPnkSt\nVmNra0vz5s1Zs2ZNmRgPHz7Mk08+iZ2dHa6urgwfPpysrCyd13v37o2joyP29vb4+vrqbchcijgj\nIXMo9E9yrn+Sc/2TnOufIXIeGRnH2LHriY6eSWzsdKKjZzJ27PqHLuQq63x3Dl2uWLECKysrdu7c\nyeLFi1mwYAGrVq0C4LfffsPDw4Np06aRlpbGxYsXtcfl5+czc+ZMlixZQmJiIp6enmRnZzN06FC2\nbdvGnj178PX1pW/fvmRkZACaK3bPPvssALt27WLZsmXMmDGDgoIC7XlzcnIICQnBwcGBPXv28Ntv\nv7Fjxw7++c9/avd55ZVXcHd3Z8+ePRw8eJAZM2ZgY2Pz4El9CDKcKoQQQpi4RYuiOXVqls62U6dm\nERExhdDQAIOd787WGS1btmT69OkANG7cmCVLlrBp0yZeeuklnJ2dMTc3x97enjp16ugcV1JSwuLF\ni2nXrp12W2Bg4B0xL+LXX39l3bp1DBo0iI0bN5KYmEhycjIeHh4ALFiwAH9/f+0xP/zwA7m5uXz/\n/ffUqFEDgK+++orAwEBOnz5No0aNOHv2LO+88w5NmzYFoFGjRhX+/I9KrsQZCZlDoX+Sc/2TnOuf\n5Fz/DJHzgoLyr9msX2+OSsUDL9HR5Z8vP9/8keJs06aNznq9evW4dOnSfY+zsLDA19dXZ9ulS5cY\nNWoUzZo1w8nJCQcHBy5dusS5c+cASExMxN3dXVvAAXTq1ElnLl1iYiJt27bVFnAAXbt2xczMjKNH\njwLwr3/9ixEjRtC7d28++ugjjh079uAf/CFJESeEEEKYOGvr4nK3h4SUoCg88BIcXP75bGxKHilO\nS0tLnXWVSqW9UeFerK2ty9z0MHToUOLj41mwYAE7d+7kwIEDeHh4UFhY+EAx3asJL8C0adM4evQo\nAwYMYMeOHbRp04Zvvvnmgd7jYUkRZyRk3or+Sc71T3Kuf5Jz/TNEzseMCcbbe5LONm/v9wkPD6oS\n56soKysrSkoqVihu376d8PBw+vTpQ4sWLahZs6bOPLoWLVpw4cIFzp8/r922e/dunaLRx8eHw4cP\nk52drd22Y8cOSktLadGihXZb48aNCQ8PZ82aNbz66qssXbr0UT5mhcmcOCGEEMLE3ZynFhExhfx8\nc2xsSggPf+qh5sNV9vluv9J1v0eGeXl5ERcXx6BBg7CysqJWrVp33bdp06Z8//33dOrUiezsbCZM\nmICVlZX29aCgIJo3b86QIUOYP38+ubm5vP3221hY3CqNBg0axLRp0xgyZAgffPABGRkZjBo1iuef\nf55GjRqRl5fH+PHjCQsLo2HDhqSnp7Nt2za6dOnywHl4GFLEGQmZt6J/knP9k5zrn+Rc/wyV89DQ\ngIcu2h7X+W6/O7W8Jrt3bvvggw8YNWoU3t7eFBYWaq/Kldc/btmyZbz22mv4+fnh7u7O9OnTuXLl\nis65f/vtN0aOHEnnzp1p2LAhn3zyCa+88op2H1tbW9avX8+4cePo1KkTNjY2DBgwgIULFwKauXiZ\nmZkMGzaMixcv4urqSr9+/fjkk08eKS8VpVKq2ZNy5eHAVU9kZByLFkVTUGCBtXUxY8YEV+o/NEII\nUV3I7zjjdK+f271eq/Jz4jIyMggKCqJp06YEBweTmZl5131LSkpo165dlXsOW2Uw1XkrkZFxvDZi\ntU6voddGrH7kJpSVwVRzXpVJzvVPcq5/knNRWap8ETd79myCgoI4fvw4vXv3Zvbs2Xfdd+HChfj4\n+BjNM8+MTWV1+75xAw4fhjVr4I3XfyU1bYHO66lpC5g4fhXXrmnughJCCCFEWVV+OLV58+bExsbi\n5uZGWloaarWapKSkMvudP3+eYcOGMWnSJObNm8eff/5Z7vnkUvPDudmd+/bmjt7ek1i4MKTM0Gdm\nJqSkQHKyZrnzv/PzoWFD8GpQTNyWseQUfVbm/SzMJmJXczaKosLL6+/9vSjz366ump5FFYlfhmyF\nEKZOfscZp4cdTq3yNzakp6fj5uYGgJubG+np6eXu9/bbbzN37lyd55mJynO37tz//vcUoqMDdIq1\nkpLbCi5PBS+XLLq1PkvDxkl4ZeyjVvJeVMeSIPo83nTgdDnv17B0MydzLMms25xk2pN8pSUpWU1I\nPuDJ9nw3kq87kXzJjsJiMxo2LFvo3VyvUwfWrtUM2d5+xe/IoXF8tRQp5IQQQhitKlHEBQUFkZaW\nVmb7rFm6RUN5d64ArFmzhjp16tCuXTuTnWsQExNj0LvI7tbt+8YNc7y8oGe3IrxUKXjlJOB8/rCm\nSDt2DGKS4Lb+OjqsrOhgloIqfyCnWKXd7E0YHcxPQEkJThcS8L2QgG/5ZyALe1LOtyY5pz0pZ1qQ\nbO7NnhIPknPrkHzNgZxCS5TSSPILyw7ZfjB11D2LOEPnvDqSnOuf5Fz/JOeislSJIm7Dhg13fe3m\nMGrdunW5ePFimeelgabx3h9//MHatWvJz88nKyuLIUOG8N1335V7zmHDhuHl5QWAk5MTvr6+2v9D\n3SwCq9r6TYZ6/5ysmw0Sb8ajed0h8w/afbIEdVoalJbe8erfezs4oG7TBpo3J8bCAjw9UYeFgZcX\n3f7zH9Lnfkbjax3Jpwa5XMTWOYc3vl8BTz5JzOrVcOUK6jp14MIFYnbtgsuXURcVwYUL7Dt/HrJ2\n0C9rRznRwTqsGURb8m9lULtH/H4XAntuplkLM8LC1LRvDwcOGCa/sq5ZP3DgQJWKpzqsHzhwoErF\nUx3Wb3rc5xfG5fafX0xMDMnJyfc9psrPiZswYQKurq5MnDiR2bNnk5mZec+bG2JjY/nkk09kTlwl\nKi6GxnUWkH7tLPnM0273JowOrONHssHMDBo1gubNdZdmzeAezRgB4iIj2RARgXl+PiU2NgSFhxMQ\nGlqx4EpL4epVOH8eLlzQXf7e5p1gx2l2lzm0AWo+Mvdir0sw8WYdOJDpRd3aJfh1tqBDF0v8/KB9\ne3B0fKB0CSGEwcjvOOP0sHPiqnwRl5GRQVhYGGfPnsXLy4uffvoJJycnUlNTGTlyJJGRkTr7x8bG\n8umnn/LHH3+Uez75gj+YoiIYNAjiN+zmP5l9WEIj8qmBDTmEk8QeH0+m//QTNG4M1taGDrdcA9v3\nJH5/3bJDtlab+LEwQ7utBDOO0Yy9dCTeqTd7LTpzMOsJ6tcuwq+TOR2622gLOwcHQ3wSIYS4NxcX\nF65du2boMMQDcnZ2JiMjo9zXjLqIq2zGWsTFGGAORX4+hIUBxUX4bHNn9o3LZfaZEhLCh1FReo3r\nQcVFRrJ4xNtkpTlqC1D7utcJXzqfAH9/OHgQ9u+HAwc0/5uQAEVFxAD+mJFEc/bSgfgaPdlr1ZVD\n2Y1wr12In58Kv4AadOiool278gs7uSv2wRjie17dSc71T3Kuf8acc6O+O1UYRm4uDBgATjWLWH75\nSXbduMwkCwtmFRdr93nf25unwsMNGGXFBISGwlJuG7J1JSh8+q0hW39/zXJTQQEcPQorV2Kel0fL\n/ftpeXA1Q7O/gxwoxpzE1BbEp/oRv74bv9h051BuYxrUzsevnYKf2p4Onc1JT49j3Bi5K1YIIcTj\nIVfiRBk3bsDTT0ND92KWXQjBIm4zeHgQN306G37++eHmrhm70lI4dUpzpe725dIlQFPYHcWHePyI\nN+vEXtse7M75AYWPy5yqU/tR/BX/pb4/gRBCCCMkw6m3kSLu3q5dgz59wLd1MZ+f6YvZpg1Qrx7E\nxWnmvYlbFAUuXtQdit2/H05rOt95ouYcW8ocZmP2NlMmfECPPvZ07Ai2tvoOXAghhLEw6menCg19\n3Dp++TL06gVdOxXzf+ef0RRwderA5s3VsoC7b85VKqhfH0JDYdIk+OUXzdW6zEyIicHStqTcwxxL\nj3N19le8E3yAWvYFdPNOY8Lwy/zxewlXrlT+5zAm0iJB/yTn+ic51z9TzbnMiROA5oLSk0/CgH4l\nzDz6PKqodZpnWm3apGkVIirO0RF69qRDcxWq/eU0Mq65j0+VWMgZTw527D7diW2ne/DZ92oG0wUP\nl1x6dCykR39XejxpwxNPVOzRYkIIIaoXGU4VnD0LvXvDsCElTDo4EH79FZydNVfgfO/2rARxP/e8\nKzYkBA4fhh07YPt2zf+mpFCMOYdowzZ6sA1/tlmowcqKHq2v06NPTXr0c6FNG7CQP7+EEKJakDlx\nt5EiTtcQ1xAyAAAgAElEQVSpU5orcGPeKuXtvYPgxx81V5I2bQI/P0OHZ/QeqJHx+fOaYu5mYbd/\nP0pJCWd44u+irgfbLNRcUHnQxfsyPXpa0ON5Nzp1s6BGDf1+LiGEEPohRdxtjLWIexw9bpKSICgI\nJr1Xyut/DYfvvoOaNWHDBujSpVLfyxgZvK9QTg7s2XPrSt2OHZCZyRVc2UE3TVFnFsBB2tKq7hV6\ndC6mx7O16R5Sk5tPp5sz/RO+XLyR0mIbzCzyGfXWk0ycPt5wn+k+DJ7zakhyrn+Sc/0z5pxLnzgB\n6DaeLSwsJikpmPmf9mDojlGaAs7ODtatkwKuqqhRA9RqzQKaNidJSdTavp1ntm/nmR1L4cRE8rBh\nT2pHtv3WgyW/9WC4WQ/c7HNxcPwfR8+fILf0VjPm2bMGA59U6UJOCCFExciVuGoiMjKOsWPXc+rU\nLO02N7f3+drvNKFrV4GNDaxdC4GBBoxSPLBLl2DnTs3Vuu3bYe9eSgqLSaAlvWjMVVaXOaSB0wBS\nMn6XmyWEEMIIyHDqbaprERcSMpno6Jllt9ORKOvD8McfEBxsgMhEpSoogPh42LEDr/fWk1K8ocwu\nZryLm/UE1D6XUPe1I3CwO42bmUtRJ4QQVZD0iTMBj9rjpqCg/JHzfGrC6tVSwJXDKPsKWVtDt24w\nfjxmjubl7uLFJrYVdKD3/rlsnRVLoE8aDeyuMNjvKEunJHPqeAmG+jvHKHNu5CTn+ic51z9TzbkU\ncdWEtXVxudtt2jWBvn31HI3Qh1FvPYmTxWCdbU4Wg3htTAiNvpnKq8NK+d5rKucUD2Lyu6DeN48t\nM7fj3/wSnnaX+YdfAsumnOH08WKDFXVCCCHuToZTq4k50z9h5swkskuWarc5mb3Mu1P8ZJK7CZsz\n/RO+WryJkmJrzC0KeO2t3mV/3ikpEBurWWJiUE6f5iSN2UIgMajZouqFlbUZap90zfDrEE+8mlga\n5gMJIUQ1I3PiblNdi7jJISFciK7HFkrw4hw25BBOErtCuvNhVNT9TyCqj3PntAUdsbEoJ09ynKaa\ngo5AYlSB2FiD2ucSgX1tUQ9tSMPGUtQJIcTjIHPiTMCjjudbFBSQSx8+BGKIJYq9hJKNeX5+pcRn\nikx1DsV9NWgAgwfD0qVw4gSq8+dptmIao0Yq/Nh0GheVuqzLV9N53+esnRlPpybXaGR3kX/6HeT7\nKcc5d6qw3NNGRsYREjIZtXo6ISGTiYyMK7NPtc25AUnO9U9yrn+mmnPpE1dNFFtbE48f05ihs73E\nxsZAEQmj4e4Or7yiWQBVaiot4uJoERPD6NjpKEkvkZjXgi37Avljn5p/zXTGwbaIQJ901H1sCRz+\nBAcS/+K1EatJTVugPe2RQ+P4aimEhgYY6pMJIYRRk+HUaiJy1Xqee6kbuThhTikA73t789TChXd/\nDJQQFZGWBnFxmuHXmBhKE5M4ig8xqLVLtupjCpRPyxzaqf0o/or/Uv8xCyGEkZA5cbeprkXc5v/d\nYNyAg/RXBWDu70+Jre29n+MpxMNKT9cUdX/PqytNOIo7fUgjssyurjZDSDi+FLcGVgYIVAghqj6Z\nE2cCHnU8P35tOoHE82GnTkyPjeXDqCgp4O7DVOdQPHZubvDii7B4MRw5gtmldOxq5JS7a16+Dc08\nc2lZM5m3uuxlxqvfcTWtSM8BV2/yPdc/ybn+mWrOpYirJuL/KsaPeGjf3tChiOqmdm06NFXwZqDO\nZm/C6Gf9G1dx5b85L9Dwr1WsXZZKo3q5tLU/zbjue/jf3ONkXim/x6EQQlR3MpxaTTSxT+P37N60\n/GocjBxp6HBENRMXGcniEW+TleZIPjWwIQf7utcJXzqfgA4dNEOvW7ZATAxFSSeJx48tBLKFQHbS\nlaYOaQT6ZhI4wBH/oY1wcJF7soQQ1YPMibtNdSzirl8Hd+ccrisOmO/dDX5+hg5JVENxkZFsiIjA\nPD+fEhubu8/JvHhRe5MEW7ZQeCKZ3XTSFnW76URLpwsEtrtO4HPO9BjSiBoO5T9iTAghjJ0Ucbcx\n1iIuJiYGtVr9UMduWZvH5NB9bLdQQ3a25vma4r4eJefi4ZSb8wsXtAUdMTHknzrPLrpoi7p9tKet\n8zkCO2QR+EItug16AtsaMlOkouR7rn+Sc/0z5pzfq26RMYlqIH5tumY+XMuWUsAJ4+PuDoMGaRbA\n5uxZ1DExqLdsYUbMUHKT09lxrRtbNgQyZUMgh0a54eeaQmDHbALDatPlJS+sbcsWdZGRcSxaFE1B\ngQXW1sWMGRMsPeuEEEZFrsRVAy/7HeepfbMYOtwcli0zdDhCVK7kZO1VOrZsIftcBtvoob1Sl6jy\noVOtMwR2ziFwoBsdX/Riw8atZZoP1687jq+WPieFnBCiSpHh1NtUxyKuqWMaq7OepNWiURAebuhw\nhHh8FAXOnNEUdX8v11Oz2Yq/tqg7qWqCyvx9soojyhwuzYeFEFWN9IkzAQ/b4+b6dUjNdqA5SdJe\n5AGZal+hquyRc65SQaNG8OqrsHw5nD+P4/G9PP3lM3z6cjz76oaSrDTEqvhiuYenJN2gsKB6/ZEn\n33P9k5zrn6nmXObEmbgDuwtpU3oQC1UptG1r6HCE0C+VCpo00SyvvQaKgsuxYzj4jeJKbtndr+XW\nxtU2h861TtOzYy49X6hNp4FPYGMnf+8KIaoeGU41cfP+dZ4z838jotlnkJRk6HCEqBIGtu9J/P66\nnGKVdps3YXQw38AXJQrb6EEsPYmlJ0fxoYNrMj39sun5rAtdB3tjW1Namggh9EPmxN2muhVxgzqf\nJGj3TIa9XAg//GDocISoEu7afHjJPAKaNdM0H/57uXHuGtvpTiw9iUHNYVrj65xCz3ZZ9BzgQrd/\neFPTSQY1hBCPh8yJMwEPO54ff7wm7dkn8+EegqnOoajK9JXzgNBQ3lo6n44hrqh7QscQV83TI55+\nWjP0OmIEfP89pKRgf/oQTy0byMdDj7HT6xXScWPatXGoNm/iwzGXqOucTxenJCb23MXaTxPJumpc\nz36V77n+Sc71z1RzbhR/PmZkZDBw4EBSUlLw8vLip59+wsnJSWefc+fOMWTIEC5duoRKpeK1115j\nzJgxBoq4arhxA85lOeLDUWjXztDhCFGlBISGlv/EiNupVPDEE5pl+HAAapw9S1BsLEGxsRD7X/JO\nnmfX9S7ExvVkblxPwsY3oLn9GXq2zqDn0/b4D2+Mc13pzyiEqHxGMZw6YcIEatWqxYQJE5gzZw7X\nrl1j9uzZOvukpaWRlpaGr68v2dnZ+Pn58fvvv9OiRQud/arTcGrclmIm9N7HLqUzXL0KLi6GDkkI\n03PhAsTFaYZfY2IoOHaG3XTSzqnbRRe8a6SjbnWFnn1r4D+8MbUa2JY5jTQfFkKUx+jnxDVv3pzY\n2Fjc3NxIS0tDrVaTdJ9J+gMGDCA8PJzevXvrbK9ORdz8iWmc/M+vfOb1iaZ3lhDi8UtLu1XUxcZS\nmHCcePy0Rd0OuuFpd5WePpfo+ZQdAcO92ZsYL82HhRDlMvo5cenp6bi5uQHg5uZGenr6PfdPTk5m\n//79dO7cWR/h6cXDjOfHb8vTPG5LhlIfiqnOoajKTCLndetCWBh89hkcOYLVpQt0/fUd3h2Tx7q2\n73GVWizLHYjX3l/4bmYKzb0LeaHf7zoFHEBq2gI+mLrisYdrEjk3MpJz/TPVnFeZOXFBQUGkpaWV\n2T5r1iyddZVKhUqluut5srOzeeGFF1i4cCE1a9as9DiNSXySHe8QD+1fNHQoQlRftWvDc89pFsAi\nI4OOW7fSMTaW8bHTKdl/CA8lhLL/+sHJIwWc3JuJt58T9/hnTwhRTVWZIm7Dhg13fe3mMGrdunW5\nePEiderUKXe/oqIinn/+eQYPHsyAAQPuer5hw4bh5eUFgJOTE76+vqjVauBWtW7s635+as5mOnCZ\nI8SYv4T6789eVeIzhnW1Wl2l4qkO6ze3VZV4Htt6//7Qv79mPTsbu5c/hWyAmL+zoNk/t/ASXTtu\nwNKiJ/5e56jbJJ62T9Zh2LgBmJlVXjw3VZn8yLqsV/K62oj+Pb/538nJydyPUcyJmzBhAq6urkyc\nOJHZs2eTmZlZ5sYGRVEYOnQorq6uzJ8//67nqi5z4rbGljK+1z7+Ku0IqalQr56hQxJC3MVdmw/X\n2MrKomskF9YjjgDiCGAr/lwxq0MPj2T8u5YQMLAe7UPrYWkll+qEMEVGPyfu3XffZcOGDTRt2pTN\nmzfz7rvvApCamkro3y0Ctm/fzvLly9myZQvt2rWjXbt2REVFGTLsSnXnX8z3E78hA7/S3Zr5OVLA\nPZQHzbl4dNU1529+OIH2dfcTQkd6oiaEjrSre4A3Vi1FdT2TJ+K+Y+jMpnwd/BPHa7QnobQF/zg7\ni7OrdvDac5dxtcnmyXpH+OCZvcR8d5a83Ir/oVpdc25IknP9M9WcV5nh1HtxcXFh48aNZbbXr1+f\nyMhIAHr06EFpaam+Q6uy4uNyUMtNDUIYhYDQUFgKGyIiMM/Pp8TGlaDw6bf62Pn7a5ZJk6CoiHoH\nDvBiXBwvxm2ErdO4dk1he1p3tv7pz3t/BnBoaC18Xc8S0C4b//4udB/khaOzUfzNLoR4AEYxnFqZ\nqstwqk+tS/xwNRjfSU/DzJmGDkcI8biUlsKRI5q2Jn8vOek3+IvO2iHYPXSkscMlAlpnEhBqj/+Q\nJ6jjbmnoyIUQFWD0feIqU3Uo4rKzwc0pn8wSeyx/XaW9K04IUQ0oCpw4oVPUFaakEo8fW/EnjgC2\n0x03uxsEtLiCf4gdAUO8aNjMRuc00nxYiKpBirjbGGsRF3PbHXv3s22rwr967Wd3sR+cPq15ZJB4\nYA+Sc1E5JOePSUoKbN2qLepKjp3gCK2II4BfcSGJ0Vhbg3+TNAJ6W1LidZ6Zc9ZJ8+HHRL7n+mfM\nOb9X3WIUc+LEg4nffB2/4r/AyQn+bqUihKjGGjbULIMHA2Cenk7brVtpGxdH68jl9Dz9IScLvIk7\nEkDckQB+5CCFlNd8eJQUcUJUIXIlzgQN6XWOgC0zGBF4GjZvNnQ4QoiqLiMDtm/XXqlruNuWs9qe\ndbfYqMYyY/BIAgY1wK+XI5YyrU6Ix87oW4yIBxN/xFrzuK327Q0dihDCGLi4QL9+MHcu/PUXFk5W\n5e7mqJziwvebGP3UaVxscniyQRIfvpRA7O/XyM/Xc8xCCCnijEVFe9zk5MCZDEdakiBF3CMy1b5C\nVZnkXP/Ky3mHJwrwZqDONm/CUNc9zkL17+y36ca5UnfePv9vbqyKZOKzx6hll0tA/RNMfjaB6JVX\nyc7W0wcwQvI91z9TzbnMiTMxBw5AS7MkrEqKpEecEOKhvPnhBBaPeJvGaR3JpwY25GBf9zpvLJ0P\noaFQUIDT7t2ExsYSGhMNO6aTnWfGzotdifs9gFm/BxCv6kDLWpcI6JhHwAt16DGgFs7Ohv5kQpgW\nmRNnYiI+zibh/eV8YfdvyMoCc3NDhySEMEJxkZG3NR+2ISg8/Fbz4TsVFsLevRAbq1m2bSM/p5jd\ndNL2qtul6koj52sE+OUQ8Gwt/J+rjZubfj+TEMZIWozcxtSLuGHBqXTfMI2R3Y5qJioLIYS+FRXB\nvn06RV1RVi77aK8t6raZBVDXPocA3ywC+rsQ8HxtGniW//xX6VknqjMp4m5jrEVcRXvctK53mW/T\nnsLvza6wePHjD8yEGXNfIWMlOdc/veS8uBgOHoSYGE1Rt3UrJZlZHKa1tqiLM1NT065U81SJpx0I\neNEN78Yq1q6N47URq02qZ518z/XPmHMufeKqidxcOHXZgVYcgfZvGjocIYTQsLAAPz/N8u9/Q0kJ\n5ocP4xsbi29sLGNiR6FkZJCU3Zy4nQFs2hnA1CmBKNY25PEz1/IidE4nPeuE0JArcSZk5054K/AI\n8QWtNUMZcmODEMIYlJZCQsKt4dfYWJTLlznDE3SkExn8WOaQeg4vcSHzR1Tlj8AKYTKkT1w1Eb89\nD7+CHWBpCS1bGjocIYSoGDMzaN0a3noLfv4Z0tNRJSTQ6PN3cLI8V+4hl7I88LTPYGjAGb5dlMXZ\ns3qOWYgqQIo4I1GRHjfxm69rmvy2bg1W5TfrFBVnqn2FqjLJuf5VyZyrVODjA6NH06GVRbk9657n\nSzbldKHb1tmsG7uODo0yaOx0hZEhZ1m5LI+0NAPFXgFVMucmzlRzLnPiTEj8AQveJF6GUYUQJqPc\nnnVumbz57oc0LSig6caNjNr2HUp+PgnXW7I5uhc/bejFm+a9qOuUTy//Inq9VIeevS1xdTX0pxGi\ncsmcOBORlweu9gVcK3HA+rP58MYbhg5JCCEqxX171uXnayYFb9qkWXbvpqQUDuDLFgLZbBbENpU/\n3m436NVLReCLtQhQm+HgYLjPJERFSYuR25hqEbdrF7wRmMi+fB/NP2Zduhg6JCGEMIzr1zU3SNws\n6hISKMKCvXRgM73YbBHCbjrSskEWgcFW9HrOie49VNjZGTpwIcqSGxtMwP3G8+N3FuJXsF0zQbhN\nG/0EZeJMdQ5FVSY51z+TzLmjIzzzDCxcCEeOQGoqlsu/pevwFkzyXM6m4p5cLnZm9pmBWH4ZwQfP\n7KGOYz4BzS8xfXw2cXFQUHDrdJGRcYSETEatnk5IyGQiI+MeKTyTzHkVZ6o5lzlxJiJ+y3U6KXuh\nRXPkz0khhLhNvXowaJBmURQ4dQqbjRtRb9qEevNCyJhGNjXYfqw7W44FMv7zUBKLm9DF5wb1Wxxk\n7Ya1XLk6T3u6I4fG8dVSpE+dMDgZTjURbT2usvTCU3Qc3By+/97Q4QghhHEoLdU8TWLTJti4EbZu\nhdxcMnFkK/6MxJF0lpc5rFP7UfwV/6UBAhbVjQynmri8PDiRZk9rDsudqUII8SDMzDT/bo4fD1FR\ncO0axMbiNHUs/bplYM35cg/bt9+Jf43OY8MG3aFXIfRJijgjca/x/EOHoJn1GWwogPbt9ReUiTPV\nORRVmeRc/yTnd7CygoAAmDEDtm/HwtGy3N3qKntx/vJjpj1/mNqOBfRTZ/H5Zwpnztz/LSTn+meq\nOZcizgTE7y7BL3+HZsXX17DBCCGECenQqLDcZsPdayYyxWI2O2604UxBfQbFvsZfE3+lS8ssmjfI\n5u3wYqKjNd1PhHhcZE6cEYuMjGPRomgOxiu4XE1krtthQtNOGDosIYQwGXGRkSwe8TZZaY63mg3X\nvU740vkE+Ptr5tGtWQORkXDpEqWo2E871lk8wzr7MA7neRPQvYQ+z9rSpw80amToTySMjfSJu42p\nFHGRkXGMHbueU6dmabd513iVhauGyh1TQghRie7bbBg0N0jEx2sKujVrYN8+ADJwZgNBrHMeRFRh\nII6uFvQZYE2fvmb07Ak2Ngb4QMKoSBF3G2Mt4mJiYlCr1dr1kJDJREfPLLNfSMgUoqI+1GNkpuvO\nnIvHT3Kuf5LzxyQ1Fdau1Vyh27ABcnIoRcUBfPnMuj3HHCZwKPsJ/P2hTz9L+vQBb+9bh98caSko\nsMDaupgxY4LlD/RHYMzf83vVLdInzkgVFJT/o8vPN9dzJEIIIcqoXx9GjNAs+fkQG4tZZCTt16zh\nH2e+Rn35a67hxIaNT7Hu0FBmvd8Dexcr+jxjiYvrVr78v9VcTF+gPZ30phPlkStxRkquxAkhhBFS\nFEhM1FyhW7MGtm+HkhJKUXGQtqxzGczHWefILl5Q5lDpTVc9SZ84EzRmTDD1647T2Va/7ljCw4MM\nFJEQQoj7UqnAxwfeeUfzfNfLl2HlSswGvUI7l7O8nzEel+ID5R56IrGEK1f0HK+o0qSIMxJ39rix\n5wbdWUtXelOD8YTQke6sw54bhgnQBJlqX6GqTHKuf5Jz/dPJubMzvPQSLF8Oly7Btm1YWJeUe1xh\nvhmNvYoJCVb45htNX2JRMab6PZcizkhFL1rET2knmEs+rXmeKPbyU9oJNkREGDo0IYQQD8PcHLp3\np4OPWbm96Z5WVnIhx5F/7nuLP2cfwcuzhH79YMUKuCF/v1dLRlPERUVF0bx5c5o0acKcOXPK3WfM\nmDE0adKEtm3bsn//fj1H+HjdeVeNxd/PecnDFlvytNvNpbNkpTHWO5mMmeRc/yTn+ne/nL/54QTa\n191PCB3piZoQOtKu1l7eCOtLjYa1GXj1c1Yfb825bGfCDk5i5axTeLiX8vzz8PPPkJurn89hTEz1\ne24URVxJSQlvvfUWUVFRHD16lJUrV5KYmKizz9q1azl58iQnTpzgq6++YvTo0QaKVj+Kra2BskVc\niTQdEkIIoxYQGspbS+fTMcQVdU/oGOJK+LcRBKxaBadPw9atMGoUDs4W/OPcR6xJbEzyjVqEJn7C\nkhkXqF9f4eWX4fff5YkRps4oirjdu3fTuHFjvLy8sLS05KWXXuJ///ufzj5//PEHQ4cOBaBz585k\nZmaSnp5uiHAfizvH84PHjGGStzd52GKH5s+u9729CQoPN0B0pslU51BUZZJz/ZOc619Fch4QGsqH\nUVFMj4nhw6ioW82FzcygRw/44gu4eBF++w2efx5nqxz+mfgO0QkeHM9vSMCJr1kw5Sr16ysMHapp\nWVdU9Hg/V1Vmqt9zoyjiLly4QIMGDbTrHh4eXLhw4b77nD9/Xm8x6ltAaCghCxey3N6N4+QxpV07\nnlq4sGwXcSGEEKbJ2hoGDIBffoH0dFiyBHr2pE7BOUbHjyDmSC2O0Bq/0z8zc+IN6tVTGDlS86Sw\n4mLNKSIj4wgJmYxaPZ2QkMlERsYZ9jOJB2IUzX5VKlWF9ruzj0pFjzMG5Y3nB4SGctT9JPuT8vhw\n9mwIDtZ/YCbMVOdQVGWSc/2TnOvfY8m5k9Ot5sJnz8LKlfD999RPSGDMtjDGAGc9e/BT8mTeHdeT\nc5dt6NAhjgMH1pOaeuvxjadOTQJMr6mwqX7PjaKIc3d359y5c9r1c+fO4eHhcc99zp8/j7u7e7nn\nGzZsGF5eXgA4OTnh6+ur/QHfvORqLOuHr53kKgrY+lSJeGRd1mVd1mW9CqxPnEhMp05w6hTqY8fg\nhx84fXYbHc4+xXjgZJvn6LmtlNSssdwSw6lTQUREbCA0NKBqfZ5qtH7zv5OTk7kvxQgUFRUpjRo1\nUs6cOaMUFBQobdu2VY4ePaqzT2RkpNKnTx9FURRl586dSufOncs9l5F85DK2bNlS7vZZ9SOUd/lI\nUfbs0W9A1cDdci4eH8m5/knO9c8gOS8uVpSNGxVl2DBFsbdXFFDa0lPRPEJCd2nSaJz+43vMjPl7\nfq+6xez+ZZ7hWVhYsHjxYkJCQvDx8WHgwIG0aNGCL7/8ki+/1DyCpG/fvjRq1IjGjRszatQoPv/8\ncwNHrR+5BRaau1NtbQ0dihBCiKrK3Bx694ZvvoG0NPjxR7Isy7/T4fQZS7p310y1uzl3TlRN8uxU\nI/dvhyXUvXGcd06NhkaNDB2OEEIII/FKq67sTvDkFKu027wJo4PZJl74x3+Zf6wPqWnmhIfDq6+C\no6MBg63G5NmpJiyvWK7ECSGEeHCN3B1YyFqdpsILWUeT0gxe+G8/tp9wY1XQUvZsL+CJJ2DcOE2b\nOlF1SBFnJG6f8Hi7vCJLTRFnZ6ffgKqBu+VcPD6Sc/2TnOtfVcl58Jgx7PB2I4q9xBBLFHvZ7u1G\n0NSp0LUrXL1KpyUjWRldi4Nhs7AuyqZTJ3juOU2/YWMa1KoqOa9sRnF3qri7vBJLuRInhBDigd3s\nKzolIgLz/HxKbGx4Kjxcs336dIiLg48+guhoGnw5mTnWHzJl8Ot86zGZf/6zFk5O8Pbb8OKLYGlp\n2M9SXcmcOGNWXEx/y0iGqb7j2ZJfwIT64gkhhKgi4uPh449h9WrN5Tdzc0oHvsyaLjOZv7ohJ07A\nW2/Ba6+Bi4umgfCiRdEUFFhgbV3MmDHBJtd3Tp/uVbfIlThjlpeneXaqVYkUcEIIIR4PPz/NrapJ\nSTBnDixfjtkPy3nmh+U888wzHPjoQ+ZvaoO3N3TrFseRI+s5e9b0GwhXBTInzkiUO57/dxFnZ1Oq\n93iqA1OdQ1GVSc71T3Kuf0ab8+bNNS1KTp6E8HCwsYE//sB3aFv+e64XR7+IIylRt4ADOHVqFhER\nGwwUtIbR5vw+pIgzZrm5mitx1lLECSGE0JOGDWHRIkhJgffeAwcH2LKFei/1pMHlQ+Uekp9vrucg\nqweZE2fMEhNp6VPKqoYTaZW8xtDRCCGEqI6uX4fPP4f58+l0uSF72FNmF4/6b3PqzHysrAwQn5GT\nPnGm6uacODuZDyeEEMJAHB01V+SSk3nCIxtvBuq87Mkgim/44OsLJjqqaTBSxBmJe82JkyLu8TDV\nORRVmeRc/yTn+meyObezo4W3W5kGwp/zB6/VnsOsD0sZMgQGD9Y8+UufTDXnUsQZs7w8crGTIk4I\nIUSVUGxtTSjZOg2EQ8mm9PQpnv3sSRKjUvDwgNatYfFiKCkxdMTGTebEGbM//8TqmRCynhqIzbrf\nDB2NEEKIai4uMpL1Y8cy69Qp7bb33dx4Ki+PgKwsqFkT5s3jaLcRvPGmihs3NNPpOnc2YNBVnPSJ\nM1ElOfkUY4F1DfkxCiGEMLy7PgWiUycYPRp+/RVeew2fp1azZflSfoh159lnoV8/TT9hFxcDfwAj\nI8OpRqK88fy8zAJsyEdlJ4/cehxMdQ5FVSY51z/Juf6Zes4DQkP5MCqK6TExfBgVpSnsateGn3+G\nH34AZ2eIikLVuhWDSr/naIKClRX4+Gja0JWWap76EBIyGbV6OiEhk4mMjHukmEw153IJx4jlZRXJ\nc1OFEEIYB5UKXn4Z1GrNM7rWrIEhQ3Dq/ysRX37J8OFujB4Nc+fGkZ29nnPn5KkP9yNz4ozY2SlL\n6DazD+fHfQrz5xs6HCGEEKJiFAW+/RbGjYOsLHB1hf/7P0qff5HWrSdz9OjMMoeEhEwhKupD/cdq\nYAeVX4cAACAASURBVNInzkTl3SjGjly5EieEEMK4qFQwfDgcPgxPPglXr0JYGGaDXqa2U/kFizz1\noSwp4oxEuXPisktkOPUxMtU5FFWZ5Fz/JOf6Jzm/jacnREdrblGtUQN+/BHrvTHl7pqX9/D9SEw1\n51LEGbG8nFIp4oQQQhg3lUpz5+rBg+DvT6/CQzjd8dQHO7ORHEvyJywMkpMNE2ZVJHPijNjmZxbw\n4Z9t2bL4KLz5pqHDEUIIIR5NaSmTfXzoeuwCETQnnxrYkEM4SWx9shd2/v9j4UJ4/XXNk75q1jR0\nwI+f9IkzUXm5ilyJE0IIYTrMzLCoW5fQY8cIZa/OS3sKM5k6Ff75T00B16wZfPQRuLjEsXhxNAUF\nFlhbFzNmTHC1uYtVhlONRHnj+bm5aIo4Ozv9B1QNmOociqpMcq5/knP9k5zfW7G1dbnbSw4fhgsX\n8PCA77/X9A3+6KM4XnxxPdHRM4mNnU509EzGjl1fpq+cqeZcijgjlpeHXIkTQghhUoLHjGGSt7fO\ntvfNzAi6dg3attX0lwO6dAEvr2gKCmbp7Hvq1CwiIjboLV5DqrTh1PT0dH755RdcXV3p378/tlJY\nVCq1Wl1mW16+6u8iro7+A6oGysu5eLwk5/onOdc/yfm9lfvorsGDCVi+HNav1zyja8wY+M9/KCgo\nv4y5sx2Jqea80q7EzZ07F3Nzc+Li4lCr1Rw5cqSyTi3uIq9AJVfihBBCmJwyj+4aPBjWroW5c8HC\nAhYtgi5dsC6+Xu7xNjYP347EmFRaERcUFMTrr7/O559/TmxsLKtXr66sUwvu0ieuwFyKuMfIVOdQ\nVGWSc/2TnOuf5PwhmZnB+PGwYwd4e8OBA4zZu5z6jqN0dqtfdyzh4UE620w155U2nHrw4EH27dtH\ncHAw7du3x8fHp7JOLe4ir8BMijghhBDVS8eOsG8fjB6N/Q8/0L3gB7LYRyItMaeIDuzBnmBDR6kX\n9+0Tl5eXV6H5bZ9++in16tVjy5Yt/PXXX1hZWTFs2DBOnz7NvHnzKi3gR2VKfeLecfyS2lmnmHDq\ndWjUyNDhCCGEEPqjKExu3ZqZCQn/3979R1dV3+kefwLnQBKqoOI45kcbDClJFpJwBQPXCwa9EjQj\nzlWraetdhiJlqYixZS4IdUmdUcHfCK2Lq2Lv8iq1ihq8lKAgRzr+ikpDXBBdNBInxkKNSoPkJOQc\n9v0DkiEkkCic7z7f73m//hn3yV5Znz52mmft/dn7SJLWa5ru17/odV2sO0pK9K9VVT4PeHKc0Hvi\n5syZo08++UQlJSWaOnWqxo4dq6SkpB7nFRcXq6WlRY8//rgk6dNPP9XmzZv15ptvnuD4OJZwR4BX\njAAAElNSkgLDh3cdnq9qva9ximqABra1+TiYOX3uxP32t7/V3//+d+3evVuvv/66Pv74Y0nSgQMH\n9MUXX3Sdd95552nKlCldxz/4wQ9UXl6up59+OgZjJ55ed+IiQW6nxpCrOxTxjMzNI3PzyPzkOfKd\ncmfoK/2jdmuH8hVNTu52nquZ93kl7qGHHlJlZaUyMzO7fT5gwAC98soramlp0dy5czVgQO998Ic/\n/OHJmRQ9UOIAAIls6ty5WlRfr7vr6yVJE/SO5p85XQtu+a8+T2ZGnztxCxYs0JIlS4758z179uiJ\nJ57QokWLTvpwseDMTlwkov8RfEX/M+kZXRl9/tAXCAMAkGC2rFun1xYv1sD339e7KfM0cNIt+n8b\nvu/3WCfN8XpLn7dT9+3bd9yfn3XWWZo+fXrMXylSVVWl3Nxc5eTkaOnSpT1+/swzz6igoEBjxozR\nBRdcoNra2pjO47twWK1KVcqgKAUOAJCwJpeW6l9fflmLJf3boK1qaHKnwPWlzxL39ddf9/lLzj33\nXG0//HRILESjUc2ZM0dVVVXasWOHVq9erbq6um7nnHPOOdqyZYtqa2t1xx136Oc//3nM5vFDj/v5\nra0KK0Upgw/6Mk8icHWHIp6RuXlkbh6Zx0BamnTKKRrz9y3atctTS0v3H7uaeZ8lbvTo0VqzZk2f\nv6gthk+CVFdXa+TIkcrKylIwGFRZWZkqKyu7nTNx4kQNHTpUklRUVKTPPvssZvPEhXD4UIlLduDW\nMAAAJyIpScrLU1ARjc1u0Xvv+T2QGX2WuJtvvlmLFi3q82u0mpubT9pQR2tqaur2YEVGRoaampqO\nef6TTz6pyy67LGbz+KHH9751ljieaYgZV79rL56RuXlkbh6Zx8jhLxmYcFaD3n23+49czbzPEjd0\n6FDdf//9mjx5slatWtXrct2uXbtiWuJ6ey/dsWzevFmrVq3qdW/OKZ0lLpV9OAAAlJcnSSoKbtU7\n7/g8iyH9+tqtyy+/XMuWLdMNN9ygu+++W9dcc43Gjx+vU089VR9++KEefPBBPfvsszEbMj09XY2N\njV3HjY2NysjI6HFebW2tZs2apaqqKp122mnH/H3l5eXKysqSJA0bNkyFhYVdLb3zvnm8HXd+1vXz\nYFBh/aNqIs1qCoV8n8/F46Oz93ueRDh+5JFHrPj/R5eOa2pqVFFRETfzJMJx52fxMo8zx5GIJGnC\nNxs15+MZ2rw5pKQk+/73vPOfGxoa1CfvW9i+fbt3+eWXe8Fg0EtKSvKSkpK8tLQ0b82aNd/m13xr\nHR0d3jnnnOPt2rXLa29v9woKCrwdO3Z0O+fTTz/1srOzvbfffvu4v+tb/keOG5s3b+7+wWuvecP0\nlfflpCt8mScR9MgcMUfm5pG5eWQeIzt3ep7kHUzP8M4+2/M++eQ/f2Rz5sfrLX2+J643e/fu1V/+\n8hclJycrNzdXgUC/LuidkPXr16uiokLRaFQzZ87U7bffrpUrV0qSZs+erRtuuEEvvfSSvv/9Q48W\nB4NBVVdX9/g9zrwn7pVXlDz9En097SdKWR/b17sAABD3olFpyBCpvV1X/tMBXfPToMrK/B7qxB2v\nt3ynEmczV0rcwdXPKfCTHyl61bVKeuF5v8cBAMB/BQVSba2W3vwf+msgU4884vdAJ+6EXvaL+HDk\nvXJJams5oMFqV1Iqj6fGytGZI/bI3DwyN4/MY+jwww0TTtne7eEGVzOnxFkq3NJx6HtTU1P9HgUA\ngPhwuMSNa39TH34otbf7PE+MUeIs0fn0SqfwvsihEseL4mLm6MwRe2RuHpmbR+YxdLjEDfnLNuXk\nSDU1hz52NXNKnKVa90UpcQAAHOnwC39VV6cJE9Tjpb+uocRZ4uj7+eFvKHGx5uoORTwjc/PI3Dwy\nj6GcHGnAAOmTT1T0Xzq69uJczTz27wZBTFDiAAA4yuDBUna2tHOn2hqf18sv71BxcUD799dr8eIB\nKi2d7PeEJxUlzhI9duL2H6TExZirOxTxjMzNI3PzyDzG8vK0budf9eDj7ygcflRvvHHo41tvXSRJ\nThU5bqdaKtzqUeIAADhaXp4eVa7q9zza7eP6+ru1fPlrPg0VG5Q4S/TYiQt7SlUrrxiJIVd3KOIZ\nmZtH5uaReYzl56tdQ476MCRJamsbaHycWKLEWSocTuJKHAAAR8vL02Dt7/VHyclRw8PEFl+7ZanH\n8x/Su3Wn6okN35emTvV7HAAA4sO+fVp3appuVanq9fuuj7OzF2rZsmnW7cQdr7fwYIOlwm0DuBIH\nAMDRTjlFpRnDpM/Wafl/+4XaBp6q5OSobrnFvgLXF26nWqLHTlw7t1Njjb0V88jcPDI3j8wNyM9X\nqb5R1b8UKxRarAULLnauwEmUOGuF2wdS4gAA6M3hr99SXZ2/c8QYJc4SR79XqLUjQImLMd7lZB6Z\nm0fm5pG5AZ0lbscOSe5mTomzVLgjeKjE8YoRAAC640oc4kmPnbgIV+Jijb0V88jcPDI3j8wN6Cxx\nH30keZ6zmVPiLBWODKLEAQDQmzPPlIYPl/btk5qa/J4mZnhPnI06OnTloFf0k6Tf6+roc1JSkt8T\nAQAQXyZPlv70J2nDBqvfp8p74lwTDiusFKUOilDgAADoxZbUVL0qKXDTTYpkZ2vq3LmaXFrq91gn\nFbdTLdHtfv7hEpcy+KBv8yQCV3co4hmZm0fm5pF57G1Zt04b3n9f/yZpcX29/vurr2rDrbdqy7p1\nfo92UlHibNRZ4liHAwCgh1cffVR3f/llt8/urq/Xa8uX+zRRbFDiLNHtHTedJS7Z8t2+OOfqe4Xi\nGZmbR+bmkXnsBdrbux0XH/6/A9vajM8SS5Q4G3ElDgCAY4oMHtzr59HkZMOTxBYlzhK97sSl8lBD\nLLG3Yh6Zm0fm5pF57E2dO1eLsrO7jkOSFmZn65JbbvFtpljg6VQbhcNqVSolDgCAXnQ+hXrH8uUa\n2Nam+tZWzbrzTueeTuU9cTZau1YpV1yi5pLrNKRqjd/TAACAGDleb+F2qoW81rDalKKU7w30exQA\nAOATSpwljtyhaGs5oEFq14AhPNkQS+ytmEfm5pG5eWRunquZU+IsFG7p4HtTAQBIcOzEWajpzv+t\n8Xf9kz6vuF96+GG/xwEAADHCTpxjwvsiXIkDACDBWVPiqqqqlJubq5ycHC1duvSY57333nsKBAJ6\n8cUXDU4Xe0fez6fEmeHqDkU8I3PzyNw8MjfP1cytKHHRaFRz5sxRVVWVduzYodWrV6uurq7X8+bP\nn69p06ZZf8v0eML7D1LiAABIcFaUuOrqao0cOVJZWVkKBoMqKytTZWVlj/OWL1+uq6++WmeeeaYP\nU8bWkd+1R4kzg+83NI/MzSNz88jcPFczt6LENTU1KTMzs+s4IyNDTU1NPc6prKzUjTfeKOnQIqCr\nwq3eoRKXmur3KAAAwCdWlLj+FLKKigotWbKk6ykO126ndtuJ6yxxXImLKVd3KOIZmZtH5uaRuXmu\nZm7Fd6emp6ersbGx67ixsVEZGRndzvnggw9UVlYmSWpubtb69esVDAY1ffr0Hr+vvLxcWVlZkqRh\nw4apsLCw61Jr57/oeDvuFAqF9P7f6pWidCklNW7m45jjk3FcU1MTV/MkwnFNTU1czZMIx53iZR6O\n4+u4858bGhrUFyveExeJRDRq1Cht2rRJaWlpOv/887V69Wrl5eX1ev6MGTN0+eWX68orr+zxMxfe\nE/dE/kN6q26YVm3IkKZO9XscAAAQI8frLVZciQsEAlqxYoVKSkoUjUY1c+ZM5eXlaeXKlZKk2bNn\n+zyhWeG2JG6nAgCQ4Ky4Ency2XolLhQKdV1yvS99mf72eYceeK9YGjfO17lcdmTmMIPMzSNz88jc\nPJsz5xsbHBM+MECpauVKHAAACYwrcRZaMPQxDW35D93+yc+lESP8HgcAAMQIV+IcE+4IsBMHAECC\no8RZ4shHj8MRSpwJR78OALFH5uaRuXlkbp6rmVPiLBSODKLEAQCQ4NiJs01Hh64eVKlrk57Xj6K/\nlxz+ejEAABIdO3EuCYcVVopSBkUpcAAAJDBKnCW67ueHw2pVqlIGH/R1nkTg6g5FPCNz88jcPDI3\nz9XMKXG26bwSxzocAAAJjZ0429TVqSD/gH73g8Ua2/CS39MAAIAYYifOJVyJAwAAosRZ48iduLBS\nlJrq6zgJwdUdinhG5uaRuXlkbp6rmVPibNPaeuhKXCpPpgIAkMjYibPN2rUacsXF2lNyvb5X9YLf\n0wAAgBhiJ84hXuvhnbhTBvo9CgAA8BElzhKd9/MP7GtXQBENTE32d6AE4OoORTwjc/PI3DwyN8/V\nzClxlgm3dPC9qQAAgJ042/x18UqN/fUV2l2xVHr4Yb/HAQAAMcROnENaWyJciQMAAJQ4W3Tezw/v\no8SZ4uoORTwjc/PI3DwyN8/VzClxlgnvP0iJAwAA7MTZZssVD2rh2iL9+29qpZtu8nscAAAQQ+zE\nOSTc6ilVrVyJAwAgwVHiLNG1Exf2uJ1qiKs7FPGMzM0jc/PI3DxXM6fEWSYcTqLEAQAAduJssyr/\nAf2p7gw9tSFdmjrV73EAAEAMsRPnkHAbV+IAAAAlzhpdO3HtAyhxhri6QxHPyNw8MjePzM1zNXNK\nnGXCBwYeKnGpqX6PAgAAfMROnGUWDv2NhrR8rkWf3CCNGOH3OAAAIIbYiXNIa0eQ26kAAIASZ4uu\nnbhIgBJniKs7FPGMzM0jc/PI3DxXM7emxFVVVSk3N1c5OTlaunRpr+eEQiGNHTtWo0ePVnFxsdkB\nDQlHuBIHAAAs2YmLRqMaNWqUNm7cqPT0dI0fP16rV69WXl5e1zl79+7VBRdcoA0bNigjI0PNzc0a\nPnx4j99l9U5cR4d+NOhl/Shpja6JrpaSkvyeCAAAxJD1O3HV1dUaOXKksrKyFAwGVVZWpsrKym7n\nPPvss7rqqquUkZEhSb0WOOuFwworRSmDohQ4AAASnBUlrqmpSZmZmV3HGRkZampq6nbOzp079dVX\nX2nKlCkaN26cnn76adNjxlQoFPrPEpds6ZVEy7i6QxHPyNw8MjePzM1zNfOA3wP0R1I/rjp1dHRo\n69at2rRpk1pbWzVx4kRNmDBBOTk5BiY0hBIHAAAOs6LEpaenq7Gxseu4sbGx67Zpp8zMTA0fPlwp\nKSlKSUnR5MmTtW3btl5LXHl5ubKysiRJw4YNU2FhYdeDEJ1tPS6P6+r0N23Xdv1NFxz+zxJX8zl2\nXFxcHFfzJMJx52fxMk+iHHeKl3k45vhkHxdb9L/nnf/c0NCgvljxYEMkEtGoUaO0adMmpaWl6fzz\nz+/xYMNHH32kOXPmaMOGDWpvb1dRUZGee+455efnd/tdVj/YsHWrcs9L1Uu5C5VX96Lf0wAAgBiz\n/sGGQCCgFStWqKSkRPn5+br22muVl5enlStXauXKlZKk3NxcTZs2TWPGjFFRUZFmzZrVo8DZLBQK\nSa2th26n8nYRI46+SoHYI3PzyNw8MjfP1cytuJ0qSZdeeqkuvfTSbp/Nnj272/G8efM0b948k2OZ\n1bkTN8SK7g0AAGLIitupJ5PVt1PXrtX3rrhIn0+doVM3PO/3NAAAIMasv52KQ7zWw1fiTrHmAioA\nAIgRSpwlQqGQOr5pV5I8BYcM8nuchODqDkU8I3PzyNw8MjfP1cwpcRYJt3QoVa18byoAAGAnzia7\n73xMBXddqT0VS6SHH/Z7HAAAEGPsxDki/E1UKQpzJQ4AAFDibBEKhRTeF6HEGeTqDkU8I3PzyNw8\nMjfP1cwpcRYJ7z9IiQMAAJLYibPKv//zA5pfOVFv/mabdNNNfo8DAABijJ04R4RbPa7EAQAASZQ4\na4RCIUqcYa7uUMQzMjePzM0jc/NczZwSZ5HWVlHiAACAJHbirPJU/v0K1f2D/s+Gs6WpU/0eBwAA\nxBg7cY4ItyVxJQ4AAEiixFkjFAop3D7g0Ndupab6PU5CcHWHIp6RuXlkbh6Zm+dq5pQ4i4QPDORK\nHAAAkMROnFUWDV2hlJbd+tUnM6URI/weBwAAxBg7cY4IdwS4EgcAACRR4qwRCoUocYa5ukMRz8jc\nPDI3j8zNczVzSpxFwtEgJQ4AAEhiJ84a6ypf18x/flanq0WZU3M0d26JSksn+z0WAACIoeP1loDh\nWfAdrFu3Rbfe9qr26AntkVT3qlRfv0iSKHIAACQobqda4NFHX1X9rmndPquvv1vLl7/m00SJwdUd\ninhG5uaRuXlkbp6rmVPiLNDe3vsF07a2gYYnAQAA8YISZ4HBgyOSint8npwcNT5LIikuLvZ7hIRD\n5uaRuXlkbp6rmVPiLHDRxGEaNvCn3T4bFvippkwY6tNEAADAb5Q4C+x7+zX9r+iLKtF4XahilWi8\n/m9krb55Z6PfoznN1R2KeEbm5pG5eWRunquZ83SqBQLt7ZqoNt2u97t9/l5bm08TAQAAv3ElzgKR\nwYN72YiTosnJpkdJKK7uUMQzMjePzM0jc/NczZwSZ4Gpc+dqUXZ2t88WZmfrkltu8WkiAADgN0qc\nBSaXluofZs7UHSUlWnzhhbqjpETTli3T5NJSv0dzmqs7FPGMzM0jc/PI3DxXM2cnzhIFEyeq+Pbb\n/R4DAADECWu+O7WqqkoVFRWKRqO64YYbNH/+/G4/b25u1nXXXafdu3crEolo3rx5Ki8v7/F7bP3u\nVAAAkHiO11usKHHRaFSjRo3Sxo0blZ6ervHjx2v16tXKy8vrOmfx4sVqb2/Xvffeq+bmZo0aNUp7\n9uxRIND9YiMlDgAA2OJ4vcWKnbjq6mqNHDlSWVlZCgaDKisrU2VlZbdzzj77bLW0tEiSWlpadMYZ\nZ/QocDZz9X5+PCNz88jcPDI3j8zNczVzK1pOU1OTMjMzu44zMjL07rvvdjtn1qxZuuiii5SWlqZ9\n+/bpD3/4g+kxAQAAjLHiSlxSUlKf59xzzz0qLCzU559/rpqaGt18883at2+fgenMcPUdN/GMzM0j\nc/PI3DwyN8/VzK24Epeenq7Gxsau48bGRmVkZHQ756233tKiRYskSdnZ2RoxYoQ+/vhjjRs3rsfv\nKy8vV1ZWliRp2LBhKiws7PoX3HnJlWOOOeaYY4455tj0cec/NzQ0qC9WPNgQiUQ0atQobdq0SWlp\naTr//PN7PNjwi1/8QkOHDtWdd96pPXv26LzzzlNtba1OP/30br/L1gcbQqFQ179omEHm5pG5eWRu\nHpmbZ3Pmx+stVlyJCwQCWrFihUpKShSNRjVz5kzl5eVp5cqVkqTZs2dr4cKFmjFjhgoKCnTw4EHd\nd999PQocAACAK6y4Ency2XolDgAAJB7rXzECAACA7ihxljhy4RFmkLl5ZG4emZtH5ua5mjklDgAA\nwELsxAEAAMQpduIAAAAcQ4mzhKv38+MZmZtH5uaRuXlkbp6rmVPiAAAALMROHAAAQJxiJw4AAMAx\nlDhLuHo/P56RuXlkbh6Zm0fm5rmaOSUOAADAQuzEAQAAxCl24gAAABxDibOEq/fz4xmZm0fm5pG5\neWRunquZU+IAAAAsxE4cAABAnGInDgAAwDGUOEu4ej8/npG5eWRuHpmbR+bmuZo5JQ4AAMBC7MQB\nAADEKXbiAAAAHEOJs4Sr9/PjGZmbR+bmkbl5ZG6eq5lT4gAAACzEThwAAECcYicOAADAMZQ4S7h6\nPz+ekbl5ZG4emZtH5ua5mjklDgAAwELsxAEAAMQpduIAAAAcQ4mzhKv38+MZmZtH5uaRuXlkbp6r\nmVPiLFFTU+P3CAmHzM0jc/PI3DwyN8/VzClxlti7d6/fIyQcMjePzM0jc/PI3DxXM6fEAQAAWIgS\nZ4mGhga/R0g4ZG4emZtH5uaRuXmuZp5wrxgpLi7WG2+84fcYAAAAfbrwwguP+WBGwpU4AAAAF3A7\nFQAAwEKUOAAAAAtR4uJMVVWVcnNzlZOTo6VLl/Z6zty5c5WTk6OCggL9+c9/Njyhe/rK/JlnnlFB\nQYHGjBmjCy64QLW1tT5M6Zb+/Pdckt577z0FAgG9+OKLBqdzU38yD4VCGjt2rEaPHq3i4mKzAzqo\nr8ybm5s1bdo0FRYWavTo0frd735nfkiH/OxnP9NZZ52lc88995jnOPf300PciEQiXnZ2trdr1y7v\nwIEDXkFBgbdjx45u56xbt8679NJLPc/zvHfeeccrKiryY1Rn9Cfzt956y9u7d6/neZ63fv16Mj9B\n/cm887wpU6Z4paWl3gsvvODDpO7oT+Zff/21l5+f7zU2Nnqe53lffPGFH6M6oz+Z33nnnd6CBQs8\nzzuU9+mnn+51dHT4Ma4TtmzZ4m3dutUbPXp0rz938e8nV+LiSHV1tUaOHKmsrCwFg0GVlZWpsrKy\n2zlr167V9ddfL0kqKirS3r17tWfPHj/GdUJ/Mp84caKGDh0q6VDmn332mR+jOqM/mUvS8uXLdfXV\nV+vMM8/0YUq39CfzZ599VldddZUyMjIkScOHD/djVGf0J/Ozzz5bLS0tkqSWlhadccYZCgQCfozr\nhEmTJum000475s9d/PtJiYsjTU1NyszM7DrOyMhQU1NTn+dQKr67/mR+pCeffFKXXXaZidGc1d//\nnldWVurGG2+UJCUlJRmd0TX9yXznzp366quvNGXKFI0bN05PP/206TGd0p/MZ82ape3btystLU0F\nBQVatmyZ6TETiot/P6n8caS/f6i8o94Kwx+47+7bZLd582atWrVKb775Zgwncl9/Mq+oqNCSJUuU\nlJQkz/N6/Hce305/Mu/o6NDWrVu1adMmtba2auLEiZowYYJycnIMTOie/mR+zz33qLCwUKFQSPX1\n9brkkku0bds2nXLKKQYmTEyu/f2kxMWR9PR0NTY2dh03NjZ23do41jmfffaZ0tPTjc3omv5kLkm1\ntbWaNWuWqqqqjnu5Hn3rT+YffPCBysrKJB1a/l6/fr2CwaCmT59udFZX9CfzzMxMDR8+XCkpKUpJ\nSdHkyZO1bds2Stx31J/M33rrLS1atEiSlJ2drREjRujjjz/WuHHjjM6aKJz8++nvSh6O1NHR4Z1z\nzjnerl27vPb29j4fbHj77bedWMz0U38y//TTT73s7Gzv7bff9mlKt/Qn8yOVl5d7a9asMTihe/qT\neV1dnXfxxRd7kUjE279/vzd69Ghv+/btPk1sv/5kftttt3mLFy/2PM/zdu/e7aWnp3tffvmlH+M6\nY9euXf16sMGVv59ciYsjgUBAK1asUElJiaLRqGbOnKm8vDytXLlSkjR79mxddtll+uMf/6iRI0dq\nyJAheuqpp3ye2m79yfyuu+7S119/3bWfFQwGVV1d7efYVutP5ji5+pN5bm6upk2bpjFjxmjAgAGa\nNWuW8vPzfZ7cXv3JfOHChZoxY4YKCgp08OBB3XfffTr99NN9ntxeP/7xj/XGG2+oublZmZmZ+vWv\nf62Ojg5J7v795Gu3AAAALMTTqQAAABaixAEAAFiIEgcAAGAhShwAAICFKHEAAAAWosQBAABYiBIH\nAABgIUocAACAhShxAAAAFqLEAcB3sH//fuXm5qqoqEiRSKTr81dffVUDBgzQY4895uN0ABIBW/kD\nQQAAAVZJREFUX7sFAN9RTU2NJkyYoNtuu0333nuv9uzZo4KCAk2cOFEvvfSS3+MBcBwlDgBOwCOP\nPKJ58+Zpw4YNuv/++7V9+3Zt27aNLzIHEHOUOAA4QaWlpdq0aZMikYhee+01TZkyxe+RACQAduIA\n4ARdd911OnDggAoKCihwAIyhxAHACdi9e7duvfVWnXfeeaqpqdGjjz7q90gAEgQlDgC+I8/zdP31\n1yslJUUbN25URUWF5s+frw8//NDv0QAkAHbiAOA7euCBB7RgwQJt3rxZkyZNUkdHhyZMmKD29na9\n//77Sk5O9ntEAA7jShwAfAdbt27Vr371Ky1cuFCTJk2SJAWDQa1evVoNDQ365S9/6fOEAFzHlTgA\nAAALcSUOAADAQpQ4AAAAC1HiAAAALESJAwAAsBAlDgAAwEKUOAAAAAtR4gAAACxEiQMAALAQJQ4A\nAMBC/x9ae7W1h5jDDAAAAABJRU5ErkJggg==\n", "text": [ "" ] } ], "prompt_number": 20 }, { "cell_type": "heading", "level": 5, "metadata": {}, "source": [ "Accuracy check" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For a closed body, the sum of all the source strengths must be zero. If not, it means the body would be adding or absorbing mass from the flow! Therfore, we should have\n", "\n", "$$\\sum_{i=1}^{N} \\sigma_i l_i = 0$$\n", "\n", "where $l_i$ is the length of the $i^{\\text{th}}$ panel.\n", "\n", "With this, we can get a measure of the accuracy of the source panel method." ] }, { "cell_type": "code", "collapsed": false, "input": [ "# calculates the accuracy\n", "accuracy = sum([panel.sigma*panel.length for panel in panels])\n", "print '--> sum of source/sink strengths:', accuracy" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "--> sum of source/sink strengths: 0.00461703117528\n" ] } ], "prompt_number": 21 }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ "Lift coefficient" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The lift is given by the Kutta-Joukowski theorem, $L = \\rho \\Gamma U_\\infty$, \n", "where $\\rho$ is the fluid density. The total circulation $\\Gamma$ is given by:\n", "\n", "$$\\Gamma = \\sum_{i=1}^N \\gamma l_i$$\n", "\n", "Finally, the lift coefficient is given by:\n", "\n", "$$C_l = \\frac{\\sum_{i=1}^N \\gamma l_i}{\\frac{1}{2}U_\\infty c}$$\n", "\n", " with $c$ the chord-length of the airoil" ] }, { "cell_type": "code", "collapsed": false, "input": [ "# calculates of the lift\n", "cl = gamma*sum(panel.length for panel in panels)/(0.5*freestream.u_inf*(x_max-x_min))\n", "print '--> Lift coefficient: Cl = %.3f' % cl" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "--> Lift coefficient: Cl = -0.013\n" ] } ], "prompt_number": 22 }, { "cell_type": "heading", "level": 5, "metadata": {}, "source": [ "Challenge task" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Based on what has been done in the previous notebook, compute and plot the streamlines and the pressure coefficient on a Cartesian grid." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---" ] }, { "cell_type": "raw", "metadata": {}, "source": [ "Please ignore the cell below. It just loads our style for the notebook." ] }, { "cell_type": "code", "collapsed": false, "input": [ "from IPython.core.display import HTML\n", "def css_styling():\n", " styles = open('../styles/custom.css', 'r').read()\n", " return HTML(styles)\n", "css_styling()" ], "language": "python", "metadata": {}, "outputs": [ { "html": [ "\n", "\n", "\n", "\n", "\n" ], "metadata": {}, "output_type": "pyout", "prompt_number": 23, "text": [ "" ] } ], "prompt_number": 23 } ], "metadata": {} } ] }