{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Working, wait for the figure, count to 30\n" ] } ], "source": [ "\"\"\" From \"COMPUTATIONAL PHYSICS\" & \"COMPUTER PROBLEMS in PHYSICS\"\n", " by RH Landau, MJ Paez, and CC Bordeianu (deceased)\n", " Copyright R Landau, Oregon State Unv, MJ Paez, Univ Antioquia, \n", " C Bordeianu, Univ Bucharest, 2020. \n", " Please respect copyright & acknowledge our work.\"\"\"\n", " \n", "# Beam.py: solves Navier-Stokes equation for the flow around a beam\n", "\n", "from numpy import * # Needed for range\n", "import pylab as p\n", "from mpl_toolkits.mplot3d import Axes3D ;\n", "from matplotlib import image;\n", "\n", "Nxmax = 70\n", "Nymax = 20; # Grid parameters\n", "u = zeros((Nxmax + 1,Nymax + 1),float) # Stream \n", "w = zeros((Nxmax + 1,Nymax + 1),float) # Vorticity\n", "V0 = 1.0 # Initial v\n", "omega = 0.1 # Relaxation param\n", "IL = 10 # Geometry \n", "H = 8\n", "T = 8 \n", "h = 1.\n", "nu = 1. # Viscosity\n", "iter = 0 # Number iterations\n", "R = V0*h/nu # Reynold number, normal units\n", "print(\"Working, wait for the figure, count to 30\")" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0\n", "1\n", "2\n", "3\n", "4\n", "5\n", "6\n", "7\n", "8\n", "9\n", "10\n", "11\n", "12\n", "13\n", "14\n", "15\n", "16\n", "17\n", "18\n", "19\n", "20\n", "21\n", "22\n", "23\n", "24\n", "25\n", "26\n", "27\n", "28\n", "29\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWYAAADnCAYAAAAtvfzfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3de7wkZX3n8c+3zzkzDBflMiJXxSRIJO6C2QmI5IJicCAqui81oGvwkkUT3dW8dCPZZL2bNa/1LgZCdASNjhIMiooC0SiKolxEREW5CDIMMA4DzAxzO6f7t39UHU7107fqPn3OqW6/79erXqer6qmqp093//rpXz31lCICMzOrjtpSV8DMzJo5MJuZVYwDs5lZxTgwm5lVjAOzmVnFTC51BczMFsKznr5H3L+pXqrsdTfuvCwiVi9wlUpzYDazsXT/pjrfv+xxpcpOHHjLygWuTl8cmM1sLAXQoLHU1RiIA7OZjaUgmI5yqYxeJK0Bng1siIgn58s+CxyRF9kbeDAijm6z7R3AFqAOzETEql7Hc2A2s7E1xBbz+cDZwCdmF0TEn84+lvRe4KEu2z89IjaWPZgDs5mNpSCoD2nIiYi4UtJh7dZJEvAi4BlDORjuLmdmY6xBlJqAlZKuLUxn9nGYPwDui4hbOqwP4HJJ15Xdr1vMZjaWAqhTusW8sUzut4PTgbVd1h8fEesl7Q9cIenmiLiy2w7dYjazsdVHi3kgkiaB/wp8tlOZiFif/90AXAwc02u/DsxmNpYCmI4oNc3DM4GbI2Jdu5WS9pC01+xj4CTgpl47dWA2s7EUBPWSUy+S1gLfBY6QtE7SK/NVp5GkMSQdJOnSfPaxwLcl/RD4PvDliPhqr+M5x2xm4ymgPqT7gETE6R2Wv6zNsvXAKfnj24Gj+j2eA7OZjaXsyr/R5MBsZmNK1NFSV2IgDsxmNpayk38OzGZmlZH1Y3ZgNjOrlIZbzGZm1eEWs5lZxQSiPqKXajgwm9nYcirDzKxCArErJpa6GgNxYDazsZRdYOJUhplZpfjkn5lZhUSIerjFbGZWKQ23mM3MqiM7+TeaIW40a21m1oNP/pmZVVDd/ZjNzKrDV/6ZmVVQw70yzMyqIxvEaDQD82jW2sysh0BMx0SpqRdJayRtkHRTYdlbJd0t6YZ8OqXDtqsl/UzSrZLOKlN3B2YzG0sRUI9aqamE84HVbZa/PyKOzqdL05WSJoCPACcDRwKnSzqy18EcmM1sTIlGyamXiLgS2DRAJY4Bbo2I2yNiF/AZ4NReGzkwm9lYCvpqMa+UdG1hOrPkYV4r6cY81bFPm/UHA3cV5tfly7ryyT8zG1t9nPzbGBGr+tz9OcA7yL4D3gG8F3hFUqZdczx67diB2czGUqAFHSg/Iu6bfSzpn4EvtSm2Dji0MH8IsL7Xvh2YzWwsBTC9gGNlSDowIu7JZ58P3NSm2DXA4ZKeANwNnAa8uNe+HZjNbExpaOMxS1oLnECWi14HvAU4QdLRZN8BdwCvysseBHw0Ik6JiBlJrwUuAyaANRHx417Hc2A2s7EUDO/Kv4g4vc3ij3Uoux44pTB/KdDSla4bB2YzG1u+g4mZWYVEyGNlmJlVSXbyz3fJNjOrEN/zz8ysUrKTf84xm5lVyqgO++nAbGZjaaGv/FtIDsxmNrZ8M1YzswqJgOmGA7OZWWVkqQwHZjOzSvGVf2ZmFeLucmZmleNUhplZ5ZS5n18VOTCb2VjKemV4rAwzs8rwBSZmZhU0qqmM0cyMjyhJz5d0l6Stkp6yiMd9iaTLF+t4S0XSWyX9y1LXY6lJOkxSSPq1bnjN9sooM1XNSAVmSb8v6TuSHpK0SdJVkn4vX/cySd9e6jr28B7gtRGxZ0T8YCEO0O5DGRGfioiTFuJ4Peqyv6S1ktbnr9lVko4trH+ZpHr+RbVV0i8kfVzSE7vs8wRJjcI2WyV9cXGekY2aRtRKTVVTvRp1IOlRZLcH/zCwL3Aw8DZgZx/7WOozAY8Het6IcYzsSXaX4P9C9ppdAHxZ0p6FMt+NiD2BRwPPBLYD10l6cpf9rs+/3Gan5yxQ/ZfEr3tLd1gixEzUSk1VU70adfZEgIhYGxH1iNgeEZdHxI2SngScCxyXt6AeBJB0vqRzJF0q6WHg6ZKWS3qPpF9Kuk/SuZJW5OX3kfQlSb+S9ED++JDZCkj6hqR35q32rZK+KGk/SZ+StFnSNZIOSyueH3Mr2V1yfyjptnx5SPqtQrnzJb0zf3yCpHWS3iBpg6R7JL28UHaFpPdKujNvjX47fx5X5kUezOt4XPprQtLT8ro+lP99WvIc35G3brdIulzSykFesIi4PSLeFxH35K/ZecAy4Ig2ZesRcVtE/CXwTeCtgxyzSNJzJf1Y0oP583pSvvzlxVa2pFslXViYvyu/+3E/+zxL0kVJ2Q9K+lD++NGSPpa/jnfn76OJfN3L8v/3+yVtavfcJR0j6bv5ce+RdLakZT3+Ba/If63cI+kNhX3V8vreJul+SRdK2rew/l8l3Zu/P66U9DuFdedL+kdJX8nfX1dJOkDSB/LPzM1axDRdL8NKZUhak38Obyos+3/5871R0sWS9u6w7R2SfiTpBknXlqn3KAXmnwN1SRdIOlnSPrMrIuKnwKvJW18RUfwHvRh4F7AX8G3gH8iC/NHAb5G1vN+cl60BHydr2T6OrPV2dlKP04CX5tv9JvDdfJt9gZ+S3da8SUTszFuFAEdFxG+WfM4HkLUkDwZeCXyk8LzfQ9YSfVp+7L8GGsAf5uv3zv8X3y3uMP8Afhn4ELAf8D6yVux+hWIvBl4O7E8WSN9Ysr5d5cFuGXBrj6L/BvzBPI/1RGAt8HrgMWR3Kf5iHsy+CfxBHqAOBKaA4/PtfoOspX9jn/tcC5yi7Jfd7K+zFwGfzje/AJghe889BTgJ+PPC7o8Fbif7n7+rzVOqA38FrASOA04E/rLHv+HpwOH5sc6S9Mx8+f8Engf8EXAQ8ADwkcJ2X8m32x+4HvhUst8XAX+X12Un2Wfg+nz+IrL31JIbco75fGB1suwK4MkR8Z/J4tPfdNn+6RFxdESsKnOwkQnMEbEZ+H2y//c/A7+SdImkx/bY9AsRcVVENMjeRP8d+KuI2BQRW4C/Jwu2RMT9EfG5iNiWr3sX2Zu36ON5y+4hsjfwbRHx7xExA/wr2YduWKaBt0fEdH4L9K3AEZJqwCuA10XE3Xlr8zsRUSat8yfALRHxyYiYiYi1wM1AMR3w8Yj4eURsBy4k+xKblzxgfRJ4W/6/62Y92ZdNJwflLcfZ6UVtyvwp8OWIuCIipsm+yFYAT4uI24EtZM/rj4DLgLsl/XY+/638/dLPPu8kC07Py8s+A9gWEVfn79GTgddHxMMRsQF4P/n7bvY5R8SH89dke3rgiLguIq7O198B/BOt783U2/Lj/Yis8XB6vvxVwN9GxLr8PfNW4AXKUygRsSYithTWHSXp0YX9XpzXZwdwMbAjIj4REXXgswz3MzAvwwrMEXElsClZdnn+uQe4GjikZcMBjVQuK28Zvwwg/xD9C/AB5t5w7dxVePwYYHeyHObsMpGlGJC0O9kHZjUw2zLdS9JE/qYDuK+wv+1t5ov50/m6v/DCA2zL978S2A24bYB9HgTcmSy7k6xVPuveNsdsIekrzLVsXxURactqttwK4IvA1RHxf0vU8WCSD0FifUT0+hA0Pc+IaEi6i7nn+U3gBLIW7DeBB8kC3XH5/CD7/DTZe/ETZL86ZlvLjydrld9TeN/VaH5vFh+3yFvr7wNWkb2HJ4Hrum2T7PNO4D8V6nOxpOKXTx14rKR7yRokLyT7vMyWWQnMfqEu5mdgYH32Y16ZpBnOy1NvZb2C7EupfVXgckkB/FOZ/Y5UYC6KiJslnU/27Q/Zk29btPB4I9kb53ci4u42Zd9Alv88NiLuzX96/wAWrDPkNrIP2awDgHUlttsI7CBLpfwwWdfp/zBrPdkHs+hxwFdLHLf5QBEn9yojaTnweeBu5l6rXp4PfKvf+iTWMxeIUBYRD83rAVnwfQ7wBLJfTQ8CLyELzGn6quw+/xV4r7LzEs/P9wVZgNwJrEy+aIt6vW7nkL0XT4+ILZJeD7ygxzaHkv0aguw1Xl+ozysi4qp0A0kvBU4lOxF7B1kq7QEW7jOwoProx7yxbJohJelvydJUbRsmwPERsV7S/sAVkm7OW+AdjUwqQ9JvKzsRdkg+fyhZ6+TqvMh9wCHdTojkP0//GXh//k9C0sGSnpUX2YsscD+Y52Jb8sVDdgPwYkkTklbT+6cp8MjzWAO8T9JB+fbH5UHwV2StnN/osPmlwBMlvVjSpKQ/BY4k6/EyVJKmyHKO24E/65AemC07IekJkj5M1pJ92zwPfyHwJ5JOzOvxBrLg+J18/TfJcrArImId2RfBarK8e6eujF33GRG/Ar5Bljb4Rf4Lj4i4B7icLGg/Ks9t/6akUq93bi9gM7A1/7X4FyW2+T+Sds9P3r2cuRbducC7JD0eQNJjJJ1aOM5O4H6yRsPf91HHSomAmUat1DQoSWcAzwZeEhFtv1wjYn3+dwNZ6ueYXvsdmcBMlhM8Fviesh4WVwM3kX04AL5O1hXtXkkbu+znTWQnn66WtBn4d+Z6CXyALGe4Md9/363IPr2OrNU221r7fB/bvhH4EVl3tE1kJzVrEbGN7KfoVXn+9anFjSLifrI30hvIPnx/DTw7Irr9zwb1tPxYJzHXS2SrpOKJveOU9VjZTBbUHgX8Xp4XHVhE/Az4b2TdKzeS/Z+fExG78vU/J8vZfyuf30x28u2qQtqqr33mPk3W2vx0svmfkZ34/AlZC/Qi4MA+ntIbydIjW8gaF51+Nhd9k+y9/jXgPRExe5HRB4FLyH5ebyF7r8/2L/8EWdrj7ryuVzPCFvICk7wx9Sbgufnnrl2ZPSTtNfuY7LNwU7uyTdt1CPJmZiNtryMOiFXnvKRU2W+c+L7ruqUyJK0l+yW3kuzX+VvIemEsJ2vgQHYO5dWSDgI+GhGn5L18Ls7XTwKfjoh2vW6ajGyO2cyslxjS5dYR0a6Dwcc6lF0PnJI/vh04qt/jOTCb2dga1UGMHJjNbCxF+NZSZmYVI+rz6HGxlBY1MK/cdyIOO3RqMQ/5iOjZTbSffS2uhfzO14j+1LPxd92NOzdGxGPms49h5ZgX26IG5sMOneL7lx26mId8RL1zF1oAGn2E2wbd9zVstS69GmvzDKwTGs0WhY2/iQNvTa9Q7Yvvkm1mVjWR5ZlHkQOzmY0t98ooYWfU+cX01qHsq1cyod7jBaknP3HSF7Db9r1+HvU69kSfWeqami9E67Z9LVk3oe7H6rsufZU2Wzrhk39mZtXjVIaZWcW4V4aZWYVEODCXsjMmuH3m0b0L5updbpLY6JHtTPO86b6mo/mpp+Vb1hde4PTY3epZxoSaM+a1JIOe5oknCutrybZpzrj3vrtvn0q370d6LLPu7u1dpAd3lzMzqxjnmM3MKiQQDffKMDOrlhFtMC92jnmKW3YeMPD29UJut5HkdVtzxBPN843Jrut39lhfnJ9pNK+bbyf2tO/xZK253/KUOs9PJWXTfbVsW2u+5VyaU+6Vs07XN5d1DtmG6Wfz29wn/8zMKmhEm8wOzGY2ttxiNjOrkAAaDQfmnnbGFL/Y2Xl41V59Dovr60lf4jTvm+aMdybrd9Sbx4XeNrOsuXw9yTnXiznmtB9zfy9+Le2XnMxPTTTnhSeTvG5x/bIkZ7ws2TZdP5XsqzWfnfR7bpnv1qe63zFARvR3po2GAIbUYpa0huyO7xsi4sn5sn3J7lZ+GHAH8KKIeKDNtmcAf5fPvjMiLuh1vNHsS2JmVkJEuamE84HVybKzgK9FxOHA1/L5JnnwfgtwLHAM8BZJ+/Q6mAOzmY2vKDn12k3ElcCmZPGpwGzr9wLgeW02fRZwRURsylvTV9Aa4FssaipjV2OCu7bNfVn06maWpjaKKYSZpDvbrnqSykhSETtmkvldU13np6eTLnEzc8eOmeT7LH1h0/n0adaaC2gi6ZLWMp+kNibn0gdTk0nqomU+SWXUOqdFsvXJsZLyk0n3u2I6oldqou9Ux6ieUreKUD8n/1ZKurYwf15EnNdjm8dGxD0AEXGPpP3blDkYuKswvy5f1pVP/pnZ+Cr/3b4xIlYtQA3afTP0rJVTGWY2ngKioVLTgO6TdCBA/ndDmzLrgOKNTg8B1vfasQOzmY0xlZwGcglwRv74DOALbcpcBpwkaZ/8pN9J+bKuFjWVMV2fYP3DnYf9THPK6XyxW9pMmlOeSXLOu5Lubsl8Y3vzvLY3f0dNJPNT03OPazPJC9njh0k6KmgkOeZkhFGS9Dn1Zc3lp5fN5X23TyX56qmke9uyND/dOV8NMJHmlJMc9ERS92L5NIfcklpP1ru7nC24Ib3FJK0FTiDLRa8j62nxbuBCSa8Efgm8MC+7Cnh1RPx5RGyS9A7gmnxXb4+I9CRiC+eYzWx8DSkwR8TpHVad2KbstcCfF+bXAGv6OZ4Ds5mNpyFeYLLYHJjNbGx5oPwSZho17n9490fm0z6G6T8xHeS6eN17vZ6s29WcmI0dzfO1bc3ll29tPvbUlub5yYeb6zKxs9Bft55UtEe/5ZYcc5pDnkpy681dqmksT3Lty+Z2WN+t+eCN5ivLqS9vXr9zeXMOeWeSY27JUU8meeEkB61CzrmW9s9Oc84t83SVljfrm8fKMDOrllH9bndgNrPxVPJy6yrq2Y9Z0hpJGyTdVFi2r6QrJN2S/+05KIeZ2eJSdvKvzFQxZVrM5wNnA58oLJsdVendks7K59/Ua0eNhti2dfkj8y3XsSf5oJYrcuqF+aQvcW1H83fM5PbuOeTlDzR/la64PxljYkvSv3fn3Ho1un8NR01d5xtTSX58MskhJznlmd06z9dXJOtWNNelvlv3fTeWJWOCJH2mG0k/6XoyjgfF+XRdmlOu9Wi+zOPz4Xy0tTWib4ueLeZ5jKpkZra0GiWnihk0x1xmVCUAJJ0JnAkwsV/nq/7MzIZqhPsxL/hYGRFxXkSsiohVE4/aY6EPZ2b2CEW5qWoGbTHfJ+nAvLXcaVSlVnXR2DLXSVctOebmWSU55uJdkmrTzesm0pzy1uZ9LXuw+b+/x4bmnPKKdc0dl2sP72iuy0xz+aKYSL7f0g666fpa83xMJjnn5c0vS3335o7NM7vP5YWn92jednr35JZbac45zUnv1ry+tc90Mq5H0gc7Cv2cW9alOeW0GdCjMRN9fGIq+NmyKhjRN8agLeYyoyqZmdkAynSXWwt8FzhC0rp8JKV3A38s6Rbgj/N5M7NKGdtURj+jKpmZVUbgS7JLqYvJLYVEZPduryS3l2saB7m2s3nd5Pbm+WWbk37Km5Kc8r3NG9TWdU+TL+SXavrWmdyjOTE8sWJ58/rd5wbEmNq9eXCM6T2bX9KZPnPQMyvSfs5JzjkZxyMKfbAbk2lOOe3P3by6ZUwRehjNz5gtpQq2hsvwJdlmNraqmKYow4HZzMaXA7OZWcU4MPemBkwVx0Fu6bfcfb7Yj3kizTFva34Flm9u3njZQzNN87Vtu7pXdqa5PPXO/ZjTXCpK+y2n67uXj+1JH+rp5rpM7Jire2378mRd8/zMjqlkfXIvxCQHPdF8aOpJP+d6knMujh0dk91zymm39Zacc2ohc8rOV4+9qva4KMN3yTaz8dVQuakHSUdIuqEwbZb0+qTMCZIeKpR586DVdirDzMbWsFrMEfEz4GgASRPA3cDFbYp+KyKePd/jOTCb2fhamFTGicBtEXHnguydxc4x12FqS2FB2m85zTG39GOe2yDNMU9ta954amsynvLD08373t68g4jm8rG9uZ9z1Avrk5yxJpJBIpJ5TSb/5okeP52SfHak+e5dc88lzT/Xkvmp6eYkcW1Xc7/n2nRz3SZ2b657fWcytkZzCrspx9xInmZM9NmPudcvSueFrR8Ll2M+DVjbYd1xkn4IrAfeGBE/HuQAbjGb2fgqH5hXSrq2MH9eRJyXFpK0DHgu8Ddt9nE98PiI2CrpFODzwOH9VTjjwGxmYyv9Fd7FxohYVaLcycD1EXFfuiIiNhceXyrpHyWtjIiNpWuRc68MM7PyTqdDGkPSAVLWF1bSMWTx9f5BDrL4/Zi3RNN80/qeY2UUcsy7mgtPJjnmya3NOeXajub5ln7JaV/iepJzni70e24ZbznJKS9blqxPv/96fB9G8o9I8sYRheeS5J+VzifPg5lGsj7NOTf3e67vlvRzTubrU8WxMpoP1To+czKf5oz7zCGP6M0pbDENMccsaXey0TRfVVj2aoCIOBd4AfAXkmaA7cBpEemHuRynMsxsPA355F9EbAP2S5adW3h8NtmNq+fNgdnMxteIXvnnwGxm48uBubdavXkMi5afGelwvvXmBWq651+SU97RnDOe2JHkWnemY1/M457lSdooGkk9G/O8H3qyv5Y0VSE/HmmuPOmPnaZha8m+lM4n//O0n3NtV5JzXjY3n/ZbbkymY4IkVe2VancO2eZB9NUro1LcYjaz8TTCgxg5MJvZ+HJg7k31YNnmuZ/eLd9mLSmB5tW1QvpBSSqjtqv5J31te3IJ9nQfw3hC76E6i6vSsrUF7h7eJVUS6fOcSP4Pade+5HmlqY60e11tunn7xmQhlTGZdDnsdWup9P+WcCrD5s2B2cysWpzKMDOrGgdmM7MKCffKKEX1YNlDhUub05xyS/e5tJtY4XLuJEesJBdaHBoTaLmsuWXfaV3ToTqL0nqnl1ynw4Cmt5rqJR1WNMkDR3H/va747DGEaLrvlu51abfCXc3PrVbMMae59XR40/R59JtD7pGTNmvhFrOZWbU4x2xmVjUOzGZmFRI4MJehmQYTmx4uv0GaPy3MK8nztvTtnUn6Kae51l455qnmf02x/2+6bZqnbckxt/SJ7pFz7jWsaHF/6f8hPVZal7R8etuqlku00+FRk5x0rZhjLt/3u9R6s3kQTmWYmVWOA7OZWdWMaGCe17XDklZL+pmkWyWdNaxKmZkNRZScKmbgFrOkCeAjZLdaWQdcI+mSiPhJx43qdXhgc8fV/Vjw/+XUso6rFjwzmuagJxdu7I30Flotw6Gm/cHNRsUIjy43n0/8McCtEXF7ROwCPgOcOpxqmZkNwa9bixk4GLirML8OODYtJOlM4EyA3Wp7zuNwZmb9GdVLsufTYm73i77luycizouIVRGxalltxTwOZ2bWH0W5qdS+pDsk/UjSDZKubbNekj6Un3O7UdLvDlrv+bSY1wGHFuYPAdZ322DzzK82XrbhnDuBlcDGeRx7IVW1blWtF7hug3Ldunv8vLZemDTF0yOi0//lZODwfDoWOIc2WYQy5hOYrwEOl/QE4G7gNODF3TaIiMcASLo2IlbN49gLpqp1q2q9wHUblOu2CBY3f3wq8InIrkC7WtLekg6MiHv63dHAqYyImAFeC1wG/BS4MCJ+POj+zMyGafbKv5KpjJWSri1MZ7bZZQCXS7quw/p2590OHqTu87rAJCIuBS6dzz7MzBZKy9ANnW0s8Qvh+IhYL2l/4ApJN0fElcXDtdlmoDb7At+crqPzlui4ZVS1blWtF7hug3LdFlLZrnIlQ2dErM//bgAuJusyXNT3ebdOliQwR0RlX/Sq1q2q9QLXbVCu28IbVq8MSXtI2mv2MXAScFNS7BLgz/LeGU8FHhokvwweK8PMxtnwTv49Frg4H0lyEvh0RHxV0qsBIuJcsrTuKcCtwDbg5YMezIHZzMbWsC7JjojbgaPaLD+38DiA1wzjeIuayqjSoEeS1kjaIOmmwrJ9JV0h6Zb87z5LVLdDJf2HpJ9K+rGk11WlfpJ2k/R9ST/M6/a2fPkTJH0vr9tnJXUebGRh6zch6QeSvlSxerVcnFCF1zOvx96SLpJ0c/6eO64qdZu3Eb0ke9ECc2HQo5OBI4HTJR25WMdv43xgdbLsLOBrEXE48LV8finMAG+IiCcBTwVek/+vqlC/ncAzIuIo4GhgdZ5P+wfg/XndHgBeuQR1A3gdWffNWVWpF2QXJxxdOPtfhdcT4IPAVyPit8lahT+tUN0Gl98lu8xUNYvZYq7UoEd5N5dNyeJTgQvyxxcAz1vUSuUi4p6IuD5/vIXsg3JwFeoXma357FQ+BfAM4KKlrJukQ4A/AT6az6sK9epiyV9PSY8C/hD4GEBE7IqIB6tQt/nqsx9zpSxmYB5a5+sF9NjZs6j53/2XuD5IOgx4CvA9KlK/PF1wA7ABuAK4DXgwv+gIlu61/QDw18BsG2i/itQL2l+cUIXX8zeAXwEfz1NAH817HVShbvMXUW6qmMUMzEPrfP3rQtKewOeA10fEcAayHoKIqEfE0WT9NI8BntSu2GLWSdKzgQ0RcV1xcZuiS/WeOz4ifpcslfcaSX+4RPVITQK/C5wTEU8BHmYU0xYduMXc29A6Xy+g+yQdCJD/3bBUFZE0RRaUPxUR/1a1+gHkP3m/QZYH31vSbC+fpXhtjweeK+kOsjTZM8ha0EtdL6DjxQlVeD3XAesi4nv5/EVkgboKdZufIV9gspgWMzA/MuhRfmb8NLIO2VVyCXBG/vgM4AtLUYk8N/ox4KcR8b7CqiWvn6THSNo7f7wCeCZZDvw/gBcsVd0i4m8i4pCIOIzsvfX1iHjJUtcLul6csOSvZ0TcC9wl6Yh80YnAT6pQt2EY1ZN/i9aPOSJmJM0OejQBrFnKQY8krQVOIBu8ZB3wFuDdwIWSXgn8EnjhElXveOClwI/yXC7A/65I/Q4ELsh72dTIBq/6kqSfAJ+R9E7gB+QnkyrgTSx9vTpdnHANS/96AvwP4FN5g+l2sgsjahWp27xUMeiWoahg4tvMbL723OfQOOrE15Uq+53P/a/rqjTMqa/8M7OxVcUTe2U4MJvZ+HJgNjOrjtkLTEaRA7OZjaeIfgbKrxQHZjMbX6MZlx2YzWx8OZVhZlYlATiVYWZWMaMZlx2YzWx8OZVhZlYxo9orY0nukm1mtuCGOLpcp9u9JWVOkPRQfvuwGyS9edCqu8VsZmMpu8BkaC3m2du9XZ+PFHidpCsi4idJuW9FxFMZLCoAAAMfSURBVLPnezAHZjMbX0MaXS6/i8vsHV22SJq93VsamIfCqQwzG1uKKDWRDf97bWE6s+M+m2/3ljouv4P8VyT9zqD1dovZzMZTf3cn2Vhm2M8et3u7Hnh8RGyVdArweeDw8hWe4xazmY2pbKyMMlMZHW73Nne0iM2zd5CPiEuBKUkrB6m5A7OZja8h3SW7y+3eimUOyMsh6Riy+Hr/INV2KsPMxlMM9dZSnW739jiAiDiX7N6SfyFpBtgOnBYD3iLKgdnMxteQustFxLfJeuB1K3M2cPYwjufAbGbjazQv/HNgNrPxpcZo3ibbgdnMxlMwtAtMFpsDs5mNJRHDvCR7UTkwm9n4cmA2M6sYB2YzswpxjtnMrHrcK8PMrFLKXW5dRQ7MZjaeAgdmM7PKGc1MhgOzmY0v92M2M6saB2YzswqJgPpo5jIcmM1sfLnFbGZWMQ7MZmYVEkDJ+/lVjQOzmY2pgHCO2cysOoKRPfnnu2Sb2fga0l2yASStlvQzSbdKOqvN+uWSPpuv/56kwwattgOzmY2vIQVmSRPAR4CTgSOB0yUdmRR7JfBARPwW8H7gHwattgOzmY2pkkG5XIv5GODWiLg9InYBnwFOTcqcClyQP74IOFFS1ztrd+LAbGbjKYBGo9zU28HAXYX5dfmytmUiYgZ4CNhvkKr75J+Zja/y/ZhXSrq2MH9eRJxXmG/X8k13XqZMKQ7MZjam+roke2NErOqyfh1waGH+EGB9hzLrJE0CjwY2la1AkVMZZjaeAiIapaYSrgEOl/QEScuA04BLkjKXAGfkj18AfD1isEsP3WI2s/E1pCv/ImJG0muBy4AJYE1E/FjS24FrI+IS4GPAJyXdStZSPm3Q4zkwm9n4GuJYGRFxKXBpsuzNhcc7gBcO41gOzGY2niLK9rioHAdmMxtfHl3OzKxKgqjXl7oSA3FgNrPx5GE/zcwqyMN+mplVRwDhFrOZWYWEB8o3M6ucUT35pwGvGDQzqzRJXwVWliy+MSJWL2R9+uHAbGZWMR7EyMysYhyYzcwqxoHZzKxiHJjNzCrGgdnMrGL+P0Vqtvhtj6oNAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWgAAADsCAYAAABZszyvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3de7xcZX3v8c9379wgREIMAgmIUinKQUGNFA89igY13sAbLXhLFcqhrR5tbSuYVlvUc/B4qvZVfampUkBFoXiQFCIQsOBBQQgqcpWbXEKQEAMkXJO99+/8sdZOZj17rnt2ZtasfN+v17z2PDPPPOuZy/7NM7/1rGcpIjAzs/IZ6ncHzMysPgdoM7OScoA2MyspB2gzs5JygDYzK6lp/e6Amdn28IbXzI7fbRhtq+71v3rmkohYsp271DEHaDOrpN9tGOXaS57bVt3hve6Yv527MykO0GZWSQGMMdbvbnTFAdrMKioYDQdoM7PSCWCE9nLQZeUAbWaVFASjA76UhQO0mVXWGA7QZmalE8CoA7SZWTl5BG1mVkIBbHEO2sysfIJwisPMrJQCRgc7PjtAm1k1ZUcSDjYHaDOrKDGK+t2JrjhAm1klBTDmFIeZWfkEsHnAl7x3gDazyhoLpzjMzEonO5LQAdrMrHQCMTrgKY7B7r2ZWRNjobYu7ZC0RNKvJd0p6eQ698+UdE5+/88kPa/b/jtAm1kljac42rm0ImkY+ArwRuBA4DhJBybVjgceiYgXAF8EPtftc3CANrNKCsSWmNbWpQ2HAndGxN0RsRn4HnB0Uudo4Mz8+nnAYkldJcEdoM2ssjoYQc+XtLrmcmLS1ELg/prymvy2unUiYgR4DHh2N/33TkIzq6QIMRptj0HXR8SiJvfXGwmnh8G0U6cjHkGbWWWNobYubVgD7FNT3htY26iOpGnArsCGbvrvAG1mlZTtJBxq69KG64D9JT1f0gzgWGBFUmcFsDS//i7gRxHdLUjtFIeZVdL4TsIpaStiRNKHgEuAYeD0iLhZ0qnA6ohYAXwT+JakO8lGzsd2u10HaDOrrNEpPNQ7IlYCK5PbPllz/WngmCnbIA7QZlZRVTiS0AHazCprrP1ZHKXkAG1mlTS+k3CQOUCbWSUFmtIcdD84QJtZJUUwZbM4+mWwe29m1lDbB6GUlgO0mVVSQCeHepeSA7SZVZZ3EpqZlVDQ/mL8ZeUAbWaV5RG0mVkJZWtxDPe7G11xgDazSgp8JKGZWWm1c77BMnOANrNKipBH0GZmZeV50GZmJeSdhGZmJZXtJHQO2syslDwP2syshHwkoZlZiY31YAQtaR5wDvA84B7gjyLikTr1RoEb8+J9EXFUq7YHe/xvZtZARHbS2HYuXToZuDwi9gcuz8v1PBURh+SXlsEZPII2s4oKxMhYT2ZxHA0ckV8/E7gC+PhUNOwRtJlV1ihq6wLMl7S65nJiB5vZIyIeBMj/PqdBvVl529dIels7DXsEbWaV1OE0u/URsajRnZIuA/asc9eyDrr03IhYK2k/4EeSboyIu5o9wAHazCpq6g71jogjG25FekjSXhHxoKS9gHUN2lib/71b0hXAS4GmAdopDjOrrLH8vIStLl1aASzNry8FLkgrSNpN0sz8+nzgcOCWVg17BG1mlRQBW3qzk/A04FxJxwP3AccASFoEnBQRJwAvAr4uaYxsYHxaRDhAm9mOqVcHqkTE74DFdW5fDZyQX/8p8OJO23aANrPKmoL0RV85B91Dkr4m6e/bqPdDSUtb1RtEku6R1HCHy45C0j9I+na/+1Fl47M42rmUlQM0IOkSSafWuf1oSb+V1PEvDUl/Iumq2tsi4qSI+HSrx0bEGyPizEbtTCVJb5Z0laRH8+f6r5Lm1Nx/haSnJW2StFHS9ZJOHt/h0aDNMyRtlvR4zeWPt9dzMGtkLIbaupRVeXvWW2cA75OUfpW+D/hORIx00thkAnof7Qp8BlhAtiNjb+DzSZ0PRcQcYC/gY8CxwMo6r1et/x0Ru9RcztkOfe8LZfy/U3Ztjp49gi6/HwDzgP82foOk3YC3AGfl5V0lnSXpYUn3Svq78X/SfJT7E0lflLSBbOGUrwGvzEePj+b1zpD0mZptHC3pl/nI9C5JS/Lbr5B0gqQXpe1IekU+73JaTTvvlPTLyTzxiDg7Ii6OiCfzBV7+lWwKUL26T0TEFcBRwCuBN09mm+MkzZT0JUlr88uXaqYiXSnpnfn1P5QUkt6Ul49s9HxbtHmrpLfU1J0mab2kl+XlwyT9NH+db5B0RE3dKyR9VtJPgCeB/eps++T8fdwk6RZJb2/xEsySdE5e/+eSDq5pa4Gk7+eft99I+h819x0q6eq8nw9K+rKkGTX3h6Q/l3RH3vanJf1e/piNks6trV9VAYzEUFuXsipvz3ooIp4CzgXeX3PzHwG3RcQNeflfyEab+wGvzut+oKb+HwB3kx3m+V7gJODqfPQ4N92mpEPJgv/fAHOBV5GthFXbr1vTdiLiOuB3wOtqqr4X+Fbnz7yuVwE3N6sQEfcBq6n5QpukZcBhwCHAwcChwN/l913JtvUNXkX22r66pnzlJNr8LnBcTd03kB1B9nNJC4GLyH5NzAP+Gvi+pN1r6r8POBGYA9xbZ9t3kb0muwL/CHxb2YELjRwN/Hu+vbOBH0iann/x/wdwA7CQbIbARyW9IX/cKPCXwHyyL8rFwJ8nbS8BXp6/Fn8LLAfeA+wDHJS8DpXkHHS1nAkcI2mnvPz+/DYkDQN/DJwSEZsi4h7gn8j+YcetjYh/iYiRPOC3cjxwekSsioixiHggIm7roK/vzfs2jyzQnN3mYxuS9DqyifafbKP6WrLA0shf5yO8RyWtb1DnPcCpEbEuIh4mC2rjr+mVFAPy/6opv5rGAbpZm2cDR0naOS+/m22v23uBlRGxMn8/VpF9Cb2ppu0zIuLm/D3ekm44Iv49Itbmjz8HuIPsC6KR6yPivLytLwCzyALqK4DdI+LUiNgcEXeT/bI5Nt/O9RFxTd6Pe4Cv17w24z4XERsj4mbgJuDSiLg7Ih4Dfkh2FFvlOUBXRERcBTwMHK3sWPlXsO2fdz4wg+Ko6V6y0c24+zvc5D60OMyziW8Db5W0C9lI//+NL9ZSS9Jza3fUNWtQ0mFkz/ddEXF7G31YCGxocv//yUf8cyNifoM6C5j4mi7Ir18N/L6kPchGw2cB+yg7CutQ4MedthkRdwK3kr12O5Olasbf433JvqDHv1QeBf6QLO8+rul7LOn9ecpq/PEHkX12GtnaXkSMAWvyvu4LLEj68glgj3w7vy/pQmU7dTcC/7POdh6quf5UnfIuzZ5LFYzPgx7kAD1IO7N64SyykfMBZCOO8Q/1emAL2T/O+NE/zwUeqHlsJG2l5dT9wO+10acJ7UTEA5KuBt5ONjr8at0HZqmIlv+Ikl5KdrjqByPi8jbq70P28/lzreq2sJbsNR1PqTw3v42IeFLS9cBHgJsiYrOknwJ/BdwVEY1G5Q3bzI2nOYaAW/KgDdn78a2I+NMm/W34nkral2yUu5gsJTWa58mb/ffvU/P4IbIdtGuBEeA3+frC9XwV+AVwXERskvRR4F1NtrPD8jzoajkLOBL4U/L0BkBEjJLlqD8raU7+z/hXZCPZRh4C9m6yM+abwAckLZY0JGmhpBd20M5ZZLnFFwPnt/Hc6pJ0EHAx8OGI+I8WdXeW9GqytQauBVZOdru57wJ/J2n3fGT8SYqv6ZXAh9iWzrgiKU+mze8Brwf+jGJaaPxXyRskDUuaJekISXu3+VxmkwXwhwEkfYBsBN3MyyW9Q9kO348CzwDXkL22GyV9XNJOeX8OkvSK/HFzgI3A4/ln5s/a7OOOJZziqJQ8n/dTsn+2FcndHwaeINtZdRXZP/fpTZr7Edko7rf1crARcS3ZTsYvAo+RBZ19O2jn/Lz++RHxRKvn1sTHgN2Bb9akQ9KdhF+WtInsy+JLwPeBJfnP8m58hizP+yuyUwH9PL9t3JVkwejHDcodt5mngq4G/ivZbJvx2+8n22n3CbIgez/ZDty2/kfydRX+KW/7IbIvzp+0eNgFZPs2HiH7JfSOiNiSDwjeSpba+Q3ZL7hvkO18hGwH5ruBTWSj9spMYZxKAYyMDbV1KStFtPolbmUl6S7gv0fEZf3ui1nZzDlgz1j01fe0VfeKxV+4vtl60P3iHPSAUjZHOMhG2GZWR5Q4fdEOB+gBpGyx7wOB901BmsGssgZ9J6ED9ACKiCP63Qezsovo6JRXpeQAbWYVJUZLvAOwHT0N0LvOG449F07v5SanzNTuSu3nt/r22yk82GMVK5vbb3pmfUTs3rpmY85Bd2DPhdP56op6M8nKb3QKw89oHxdnGdb2S1kPb8fgbzuexfvdXm+9k7Z1eFbvUnKKw8yqKbI89CAb7ASNmVkTvTirt6RjJN0saUzZiWIb1Vsi6deS7pR0cjttewTdQCcpjVYpi7EyfQ92OKLoJCXS6jVzCsR6KehZDvom4B1kqwrWla+I+RWyZYLXANdJWtHqzN4O0GZWUWJ0rCdn9b4VQE1PMMShwJ350rFI+h7Z0gIO0Ga2Y+pgBD1f0uqa8vKIWD6FXVlIcbnaNWQn+WjKAdrMKimiowC9vtlaHJIuA/asc9eyiLigjfbrdaRlzs8BOtcqf9osz5zmmLtpq157zQzR2bS54eRnWMu8cAdp41b56vR1cU7atrepmmYXEUd22cQaatb/Ztva302VaO+VmdnUimjv0gPXAftLen6+tvuxTFzSeAIHaDOrrAi1demGpLdLWkN2At+LJF2S375A0sqsHzFCdrKJS8hOu3Zufr7IppziMLNKGj8n4XbfTsT51DmrUUSspeakwxGxkg7PQrTDBuhO88S1eeH0sVui+DKmj221rbEpPPR7KMkDp3nesSRnPdoih90qZ12bA0+fd6c56Qnbdo7autHZTsJS2mEDtJntAAb8O94B2swqyyNoM7OSGvTFknaYAN1pznkLw8VyTZ55czS+L2ureY56Yt8mn4MeTnLIad43nSc9rOIndrpGkvY6y1nX5qgn5IyT4vZc6tQsFQHhBfvNzMrJI2gzs7JygC6nTlMaT8f0puUnx2Zuu2+sed102lzal06n1TWbOtdqWl16/3SNFsozNNz0/uFo/nioSZFMSGEkz7PDlIcPDbfudH8QSr9VNkCbmQ36d7oDtJlVkw9UMTMrMQfowdBqGt2msZ0K5YdH5hTKj43svPX6k2MzCveNjKU55+Y55nRqXGpIaR65WK59/PShYk44zRF3Wp45tKVQnpFMw0ufW22eeAZJfjp5rHPS1nMD/hHZYQK0me2AHKDNzEoocIrDzKysfKBKSbRc0jPJfz5RM68Z4M6n9yiU73tqXqG8aWRb/c2jxZctnWuc5oynDRXvnzahfvPDsdP6tXnnoeQ3XJpDTrc9K7l/ZpIn3pIcxp62N0vFckdHqU9xTtqsJQdoM7Ny0phTHGZm5RMM/Ai65Q9USadLWifppprb5klaJemO/O9u27ebZmadUraTsJ1LN1uRjpF0s6QxSYua1LtH0o2SfilpdTtttzOCPgP4MnBWzW0nA5dHxGmSTs7LH29ng73Saq2Nh0eeVShfsW7/Qnn947ML5dojkmZMK+ZSZ04vlmel9w8XyzOG0/UsitK88rRkrnNtXnlifro4R3vGUHHbOw0X3/KZQ81z0BPWEVGT7/Tkromnx0pPidX8dWiltm+eE2119eZjcRPwDuDrbdR9TUSsb7fhliPoiPgxsCG5+WjgzPz6mcDb2t2gmVnPRJuXbjYRcWtE/Lq7Vuqb7GrWe0TEgwD53+c0qijpREmrJa1+dEN3IyYzs7YFMKb2Lr3r0aWSrpd0YjsP2O47CSNiObAc4IAXz/LvUDPrGbUfceYneeHleezK2pEuA/as87hlEXFBm9s4PCLWSnoOsErSbXmGoqHJBuiHJO0VEQ9K2gtYN8l2Jq3VvOd0rY103vMdTxXnPa+9ZkGx/VnJ9uZv3lbY9enCfTvNSNavSHLGs5IcdDrv+bHNxXVAUhPW4qh5fJrPTnPOM4aKb/HmsWJ5p+Fi38eGk5xz+joPF5/7UM160ena0a1Ox5WePsvzom3KtR+g10dEwx18EXFk112JWJv/XSfpfOBQoGmAnmyKYwWwNL++FGj3G8TMbIcjabakOePXgdeT7Vxsqp1pdt8FrgYOkLRG0vHAacDrJN0BvC4vm5mViqK9S1fbkN4uaQ3wSuAiSZfkty+QtDKvtgdwlaQbgGuBiyLi4lZtt0xxRMRxDe5a3Fbvzcz6pQeLJUXE+cD5dW5fC7wpv343cHCnbVfmSMJ03nNanrD2xqbdC+Xn/f3VhfK9p76yUN5p7lNbrz9/XnHW4R6zNhXKs6c9Uyinay4/umXnQvmmtXsVys/e9YlCeWaSZx4b2vahS5/n5uQcg2mOemTC69TZBzido107/7jVWtNpX6d6XrRZQUCLpddLrzIB2sws1W36ot8coM2suhygzcxKygG6nNL1n8eSXOtLdn2gUH757fcUynOHb5z0ttNca9qX1Fvn/aJQvnLjCwvlB56eWyg/vmVbPn3LaLJ2RpqT7jCtO3Et6uZ55VlD27afruORlqcn60EXV0cxm1pTMUOj3yoboM3MengY93bhAG1mleURtJlZWTlAl0Oa5221HvSdTxTnQX/36sMK5YP/y73F8tw1W6/vNf3Rwn1zh58slHceKs6Dnp2U5w49VSjPG9pcKG+YvUuhnK7F8aB23Xp945bioiHPjBTf0khy7yNjxddlJFnfectYMW/8TIu1O2pf57HobOWACe9ZMmm1k9bStVm8PrThHLSZWYk5QJuZlZQDtJlZOTnFMSDSnOSu04vrGu+292OFcrpm8yM162fMTObzptLHPovituYMFfO4C6bNTMqPFMoPz5hTKD8xuq3+5iRnnM73TnPOqTRvnD5+Yv2pm7Y0NIULJTjnbHUN+MdihwnQZraD8U5CM7MSc4A2MyspB+hySPOZ6fnrZiZ5391nFNdw3m+33xXKz56ZrMlcc66/9Dx9oxPW/UjWPU76Mj351MxUcY72s4bSnHWxPHt427zqx4dnJNsu9m2zim9xup7zjOR8idOG0rU3kr4na3HUvs7p8+w0x+w8sk0l4RSHmVk5BQz6eYYdoM2sujyCLqf053V6uPVu04spjH13Lp7GKj2keVZNiiRNl6Q/+9Of+q2MRvM0Qnro+C41KY6nphVTHOlyoVvGmk8JnJ6kNGYPFw87T0/ftXNyf+1rMUvF12VG8jzScpqGMptyPQjQkj4PvBXYDNwFfCAiHq1Tbwnwz8Aw8I2IaHmy7c4WTzAzGyC9OKs3sAo4KCJeAtwOnDKhH9Iw8BXgjcCBwHGSDmzVsAO0mVVXtHnpZhMRl0bE+E/Va4C961Q7FLgzIu6OiM3A94CjW7XtAG1m1ZTvJGznAsyXtLrmcuIkt/pB4Id1bl8I3F9TXpPf1lRlctBpPnMGyamZkvzovOFiDnp4ZvFrtFleecJpn7Q5KRe3ldqSTNN7PIp53tHkbZmV5LznTSv2vdbssWJOutWh3tOGiq9buq2dk6VQJ5a39T197ITXMJ0KmQxdWk2VTHlanrXU/kdkfUQsanSnpMuAPevctSwiLsjrLANGgO/Ua2IyvatMgDYzS03VPOiIOLLpdqSlwFuAxRFRb6trgH1qynsDa1tt1ykOM6uuHuSg89kZHweOiognG1S7Dthf0vMlzQCOBVa0atsB2syqqd3g3P0o+8vAHGCVpF9K+hqApAWSVgLkOxE/BFwC3AqcGxE3t2p4YFMcaf4xPeVRms+cleROn63m86RTtTnodNsT8t9J7jW1aWx6cksxd7slikuIpn3bfdrGrdfTOdLpY9PD0FPDNJ+DPT1ZWjV9brX107rp6zKd5q+Lc842lUT9xO9Ui4gXNLh9LfCmmvJKYGUnbQ9sgDYza2XQj4VygDaz6hrwH11d5aAlLZH0a0l3Sjp5qjplZjYlepOD3m4mPYKuOXTxdWRTSK6TtCIibpmqznViYl44WQcizX8mc5e35+7SdPnRTTGrWB4rllOzla6PsblBTTPbqgJnVOkmLE3q0EUzs54Z8BF0NwG6rUMXJZ04fvjkoxua78U3M5tKPVosabvpJkC3dehiRCyPiEURsWjuvOE6DzEz2z46WIujlLqZxdHxoYu33/TM+sX73X4vMB9Y38W2t6ey9q2s/QL3bbLct+b27erRJU9ftKObAL310EXgAbJDF9/d7AERsTuApNXNFibpp7L2raz9Avdtsty3HthRA3REjEgaP3RxGDi9nUMXzcx6YYc/aexkDl00M+uZHTlAd2F5n7bbjrL2raz9Avdtsty37SlAY4MdoVV/6VIzs8E2e/4+ceBRf9lW3dX/9rHry5hz91ocZlZdAz7+dIA2s8oa9J2EPV2wv0yLK0k6XdI6STfV3DZP0ipJd+R/d+tT3/aR9J+SbpV0s6SPlKV/kmZJulbSDXnf/jG//fmSfpb37Zz8rBE9J2lY0i8kXViyft0j6cZ8QffV+W19fz/zfsyVdJ6k2/LP3CvL0reu7cCHenekZnGlNwIHAsdJOrBX26/jDGBJctvJwOURsT9weV7uhxHgYxHxIuAw4C/y16oM/XsGeG1EHAwcAiyRdBjwOeCLed8eAY7vQ98APkJ2xopxZekXwGsi4pCaXGcZ3k+AfwYujogXAgeTvX5l6dvktXmYd5lH2b0cQZdqcaWI+DGwIbn5aODM/PqZwNt62qlcRDwYET/Pr28i+4dZWIb+RebxvDg9vwTwWuC8fvZN0t7Am4Fv5GWVoV9N9P39lPQs4FXANwEiYnNEPFqGvnVLDP6h3r0M0G0trtRne0TEg5AFSeA5fe4Pkp4HvBT4GSXpX55G+CWwDlgF3AU8mp93Dfr33n4J+FvYeh6vZ5ekX5B9iV0q6XpJJ+a3leH93A94GPi3PDX0DUmzS9K37kW0dympXgbothZXsm0k7QJ8H/hoRGxsVb9XImI0Ig4hW3/lUOBF9ar1sk+S3gKsi4jra2+uU7Vfn7nDI+JlZCm+v5D0qj71IzUNeBnw1Yh4KfAEg5jOaKAXKQ5Jn8/z97+SdL6kuQ3qTdgP0UovA3THiyv1wUOS9gLI/67rV0ckTScLzt+JiP9btv4B5D+FryDLk8+VND4rqB/v7eHAUZLuIUufvZZsRN3vfgFbTyBKRKwDzif7YivD+7kGWBMRP8vL55EF7DL0rTu9O6v3KuCgiHgJcDtwSpO66X6IpnoZoLcurpTvST8WWNHD7bdjBbA0v74UuKAfnchzp98Ebo2IL9Tc1ff+Sdp9fIQgaSfgSLIc+X8C7+pX3yLilIjYOyKeR/bZ+lFEvKff/QKQNFvSnPHrwOuBmyjB+xkRvwXul3RAftNi4JYy9G0q9CIHHRGX1qTRriEbCEyJns2DLtviSpK+CxwBzJe0BvgUcBpwrqTjgfuAY/rUvcOB9wE35rlegE+UpH97AWfms3KGgHMj4kJJtwDfk/QZ4BfkO51K4OP0v197AOdn37tMA86OiIslXUf/30+ADwPfyQdOdwMfIH9vS9C37rQ/Op6fpB2WR8RkDnf/IHBOk95cKimAr7fTvg/1NrNK2mW3feLgxR9pq+5Pv/83TQ/1lnQZsGedu5ZFxAV5nWXAIuAdUSewSloQEWslPYcsLfLhfDZZQz6S0Mwqa6rmOEfEkU23Iy0F3gIsrhec8za27oeQNL4fommA7umRhGZmPdWDnYSSlpCl0o6KiCcb1Gm0H6IpB2gzq6TxBft7cCThl4E5wKp8Ct3XIEtpSBpfL38P4CpJNwDXAhdFxMWtGnaKw8yqqUcHoUTECxrcvhZ4U379brLD6DviAG1mlVXmw7jb4QBtZpVV5oWQ2uEAbWbVFMCAn/LKAdrMqmuw47MDtJlVl1McZmZlNeBHSjtAm1k1hWdxmJmVUnagikfQZmbl5BG0mVk5eQRtZlZGU3O2lL5ygDazigrkA1XMzErKKQ4zsxLyNDszsxLzCNrMrKQGOz47QJtZdXmanZlZGQUw6gBtZlY6IgZ+BO2TxppZdY2fl7DVpQuSPi3pV/kJYy+VtKBBvaWS7sgvS9tp2wHazKqrBwEa+HxEvCQiDgEuBD6ZVpA0D/gU8AfAocCnJO3WqmEHaDOrpiBbLKmdSzebidhYU5xN/bkjbwBWRcSGiHgEWAUsadW2c9BmVlkd5KDnS1pdU14eEcvb3o70WeD9wGPAa+pUWQjcX1Nek9/WlAO0mVVUwFjbw+P1EbGo0Z2SLgP2rHPXsoi4ICKWAcsknQJ8iCydUWiifgebc4A2s2oKpuxIwog4ss2qZwMXMTFArwGOqCnvDVzRqjHnoM2sunqQg5a0f03xKOC2OtUuAV4vabd85+Dr89ua8gjazCqrR/OgT5N0AFmovxc4CUDSIuCkiDghIjZI+jRwXf6YUyNiQ6uGHaDNrLp6EKAj4p0Nbl8NnFBTPh04vZO2HaDNrJoiYHSw1xt1gDaz6hrwQ70doM2suhygzcxKKACfk9DMrIwCwjloM7NycorDzKyEAs/iMDMrLY+gzczKaErWeu4rB2gzq6agk9XsSskB2syqyyNoM7OScoA2MyuhCGJ0tN+96IoDtJlVl48kNDMrKac4zMxKKDo6J2EpOUCbWXV5BG1mVkbeSWhmVk5ebtTMrMR6sNxofjLYo8lOGrsO+JOIWFun3ihwY168LyKOatX20FR21MysLAKIsWjr0qXPR8RLIuIQ4ELgkw3qPRURh+SXlsEZPII2s6qK3izYHxEba4qzyb4bpoQDtJlVVgej4/mSVteUl0fE8nYfLOmzwPuBx4DXNKg2K9/GCHBaRPygZbsx4NNQzMzqkXQxML/N6usjYkmTti4D9qxz17KIuKCm3inArIj4VJ02FkTEWkn7AT8CFkfEXU2fgwO0mdnUkLQvcFFEHNSi3hnAhRFxXrN63kloZtYFSfvXFI8CbqtTZzdJM/Pr84HDgVtate0ctJlZd06TdADZNLt7gZMAJC0CToqIE4AXAV+XNEY2MD4tIloGaKc4zMxKykMQsKgAAAAsSURBVCkOM7OScoA2MyspB2gzs5JygDYzKykHaDOzknKANjMrKQdoM7OS+v8hqQu1MWoyhwAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "finished\n" ] } ], "source": [ "\n", "def borders(): # Initialize stream,vorticity, sets BC\n", " for i in range(0, Nxmax+1): # Initialize stream function\n", " for j in range(0, Nymax+1 ): # Init vorticity\n", " w[i,j] = 0.\n", " u[i,j] = j*V0\n", " for i in range(0,Nxmax+1 ): # Fluid surface\n", " u[i,Nymax] = u[i,Nymax-1] + V0*h\n", " w[i,Nymax-1] = 0. \n", " for j in range(0,Nymax+1 ):\n", " u[1,j] = u[0,j]\n", " w[0,j] = 0. # Inlet\n", " for i in range(0,Nxmax+1 ): # Centerline\n", " if i <= IL and i>=IL + T:\n", " u[i,0] = 0.\n", " w[i,0] = 0. \n", " for j in range(1, Nymax ): # Outlet\n", " w[Nxmax,j] = w[Nxmax-1,j] \n", " u[Nxmax,j] = u[Nxmax-1,j] # Borders\n", "\n", "def beam(): # BC for the beam\n", " for j in range (0, H+1): # Beam sides\n", " w[IL,j]=-2*u[IL-1,j]/(h*h) # Front side\n", " w[IL + T,j]=-2*u[IL + T + 1,j]/(h*h) # Back side\n", " for i in range(IL,IL + T+1):\n", " w[i,H-1] = -2*u[i,H]/(h*h); # Top\n", " for i in range(IL, IL + T+1 ):\n", " for j in range(0,H+1):\n", " u[IL,j] = 0. # Front\n", " u[IL + T,j] = 0. # Back\n", " u[i,H] = 0; # Top\n", " \n", "def relax(): # Method to relax stream\n", " beam() # Reset conditions at beam\n", " for i in range(1, Nxmax): # Relax stream function\n", " for j in range (1, Nymax):\n", " r1 = omega*((u[i+1,j]+u[i-1,j]+u[i,j+1]+u[i,j-1]+h*h*w[i,j])/4 -u[i,j]) \n", " u[i,j] += r1 \n", " for i in range(1, Nxmax): # Relax vorticity\n", " for j in range(1,Nymax):\n", " a1 = w[i + 1,j] + w[i-1,j] + w[i,j + 1]+w[i,j-1]\n", " a2 = (u[i,j + 1]-u[i,j-1])*(w[i + 1,j]-w[i-1,j])\n", " a3 = (u[i + 1,j]-u[i-1,j])*(w[i,j+1]-w[i,j-1])\n", " r2 = omega*((a1-(R/4.)*(a2-a3))/4.0 -w[i,j])\n", " w[i,j] += r2\n", "\n", "m = 0\n", "borders()\n", "while (iter <= 300):\n", " iter += 1\n", " if iter%10 == 0:\n", " print(m)\n", " m += 1\n", " relax()\n", "for i in range (0, Nxmax+1):\n", " \n", " for j in range(0,Nymax+1 ): \n", " u[i,j] = u[i,j]/(V0*h) #stream in V0h units\n", "#u.resize((70,70));\n", "#w.resize((70,70));\n", "x=list(range(0,Nxmax-1)) #to plot lines in x axis\n", "y=list(range(0,Nymax-1))\n", "#x=range(0,69) #to plot lines in x axis\n", "#y=range(0,69)\n", "X,Y=p.meshgrid(x,y) #grid for position and time\n", "\n", "def functz(u): #returns stream flow to plot\n", " z = u[X,Y] #for several iterations\n", " return z\n", "\n", "def functz1(w): #returns stream flow to plot\n", " z1 = w[X,Y] #for several iterations\n", " return z1\n", " \n", "Z = functz(u)\n", "Z1 = functz1(w)\n", "fig1 = p.figure()\n", "p.title('Stream function - 2D Flow over a beam')\n", "p.imshow(Z, origin='lower');\n", "p.colorbar();\n", "fig2=p.figure()\n", "p.title('Vorticity - 2D Flow over a beam')\n", "p.imshow(Z1, origin='lower');\n", "p.colorbar();\n", "p.show() # Shows the figure, close Python shell to\n", " # Finish watching the figure \n", "print(\"finished\") " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true, "deletable": true, "editable": true }, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 2", "language": "python", "name": "python2" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.13" } }, "nbformat": 4, "nbformat_minor": 2 }