Skip to content
Permalink
master
Switch branches/tags

Name already in use

A tag already exists with the provided branch name. Many Git commands accept both tag and branch names, so creating this branch may cause unexpected behavior. Are you sure you want to create this branch?

ME3263-Lab_04/ME3263-Lab_04.ipynb

Go to file
 
 
Cannot retrieve contributors at this time
444 lines (444 sloc) 42.5 KB
Raw Blame
Open in GitHub Desktop
{
"cells": [
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"\n",
"import matplotlib.pyplot as plt\n",
"import pretty_plots # script to set up LaTex and increase line-width and font size\n",
"from scipy import linalg\n",
"\n",
"pretty_plots.setdefaults()\n",
"pi=np.pi"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# ME 3263 Introduction to Sensors and Data Analysis (Fall 2018)\n",
"\n",
"## Lab #4 Predicting Natural Frequencies with the Finite Element Method\n",
"\n",
"\n",
"### What is the Finite Element Method?\n",
"\n",
"The Euler-Lagrange dynamic beam equation is an example of a partial differential\n",
"equation (PDE). These equations are common in many engineering applications e.g.\n",
"solid mechanics, electromagnetics, fluid mechanics, and quantum mechanics. The\n",
"finite element method solves PDEs. The FEM process involves two steps to create\n",
"matrices for a computer algorithm solution. First, the PDE is integrated from\n",
"the strong form to the weak form. Second, an approximation of the variable\n",
"\"shapes\" within each \"element\" is created to convert the integrals and\n",
"derivatives into matrices\n",
"[(1)](http://bcs.wiley.com/he-bcs/Books?action=index&bcsId=3625&itemId=0470035803).\n",
"For elements with nodes only at vertices, such as cubes (hexahedrons) or\n",
"pyramids (tetrahedrals), the \"shape\" function is linear for displacement. The\n",
"simplest FEM divides a continuous structure into finite, interacting elements.\n",
"As a simple example, consider a rectangular beam with a force applied to stretch\n",
"it axially, as seen in Figure 1. Hooke's law relates the force applied to strain\n",
"and deflection as such \n",
"\n",
"$F^{e} = A^{e}E^{e}\\epsilon^{e} =\n",
"A^{e}E^{e}\\frac{\\delta^{e}}{l^{e}}=k^{e}\\delta^{e}$\n",
"(1)\n",
"\n",
"where $k^{e}$ is the stiffness (force/displacement) of each element, $E^{e}$ is\n",
"the Young's modulus of each element, $\\epsilon^{e}$ is the strain in each\n",
"element, $\\delta^{e}$ is the change in length of each element, and $l^{e}$ is\n",
"the original length of each element. Considering the displacement of each node,\n",
"$u_{i}$, the force on the element 1 is formed as such\n",
"\n",
"$\\left[\\begin{array}{c} F_{1} \\\\ F_{2}\\end{array}\\right] =\n",
"\\left[\\begin{array}{c c} k^{e} & -k^{e}\\\\ -k^{e} & k^{e}\\end{array}\\right]\n",
"\\left[\\begin{array}{c} u_{1} \\\\ u_{2}\\end{array}\\right]$ (2)\n",
"\n",
"where $F_1$ and $F_2$ are the forces on nodes 1 and 2, respectively, and $u_1$\n",
"and $u_2$ are the displacements of nodes 1 and 2, respectively. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"![Figure 1: Diagram of axial loads on a beam made of one finite element.](./figure_01.png)\n",
"\n",
"*Figure 1: Diagram of axial loads on a beam made of one finite element.*"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Predicting Natural Frequencies\n",
"\n",
"In this lab you will use Ansys to predict natural frequencies of the rectangular\n",
"beam and turbine blade. In finite element analysis (FEA), this is called a modal\n",
"analysis. In any FEA simulation, there are three steps: pre-process, solve, and\n",
"post-process. Here are the steps in a FEA modal analysis [(2)](https://ansyshelp.ansys.com/account/secured?returnurl=/Views/Secured/corp/main_page.html):\n",
"\n",
"1. Pre-process\n",
" \n",
" a. Define material constants\n",
"\n",
" b. Create geometric model\n",
"\n",
" c. Define boundary conditions\n",
"\n",
" d. Mesh the model (turn the continuous volume into discrete nodes and elements)\n",
"\n",
"2. Solve\n",
"\n",
" a. Assemble stiffness matrix, $K$\n",
"\n",
" b. Assemble mass matrix, $M$\n",
"\n",
" c. solve the eigenvalue problem, $Ku=\\lambda M u$\n",
"\n",
"3. Post-process\n",
" \n",
" a. Determine natural frequencies\n",
"\n",
" a. Plot mode shapes\n",
"\n",
" b. Determine maximum displacement/acceleration for each mode\n",
"\n",
" c. Determine maximum stress and strain for each mode"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Two Degree of Freedom Modal Analysis Example\n",
"\n",
"The solution for a modal analysis is accomplished by solving for eigenvalues of\n",
"the mass and stiffness matrices, $\\bar{\\bar{M}}$ and $\\bar{\\bar{K}}$, respectively [(3)](https://www.scribd.com/document/270286438/Engineering-Vibration-Inman-D-J). The output is a\n",
"number of eigenvalues, or natural frequencies, and their corresponding\n",
"eigenvectors, or mode shapes. As a simple example, consider the lumped mass\n",
"solution of two masses connected by three springs as seen in Fig. 2. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"![Figure 2: Two masses connected to 3 springs](./spring_mass.png)\n",
"\n",
"*Figure 2: Two masses connected to 3 springs*"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In the 2-mass system, we have 2 degrees of freedom, so there are 2 differential\n",
"equations that describe the motion of masses as such\n",
"\n",
"$m_1 \\ddot{x}_1 = -k_1x_1-k_2(x_1-x_2)$ (3a)\n",
"\n",
"$m_2 \\ddot{x}_2 = -k_3x_2-k_2(x_2-x_1)$ (3b)\n",
"\n",
"where masses 1 and 2 have mass $m_1$ and $m_2$, respectively, and springs 1, 2,\n",
"and 3 have stiffness $k_1$, $k_2$, and $k_3$, respectively. The differential\n",
"equations relate acceleration of mass 1 and mass 2, $\\ddot{x}_1$\n",
"and$\\ddot{x}_2$, to displacement, $x_1$ and $x_2$. The mass and stiffness\n",
"matrices for this problem becomes\n",
"\n",
"$\\bar{\\bar{M}}\\frac{d^2\\bar{x}}{dt^2}=-\\bar{\\bar{K}}\\bar{x}$ (4a)\n",
"\n",
"$\\left[\\begin{array}{cc}\n",
"m_1 & 0 \\\\\n",
"0 & m_2 \\end{array}\\right]\n",
"\\frac{d^2}{dt^2}\\left[\\begin{array}{c}\n",
"x_1 \\\\\n",
"x_2 \\end{array}\\right]=\n",
"-\\left[\\begin{array}{cc}\n",
"k_1+k_2 & -k_2 \\\\\n",
"-k_2 & k_2+k_3 \\end{array}\\right]\n",
"\\left[\\begin{array}{c}\n",
"x_1 \\\\\n",
"x_2 \\end{array}\\right]$ (4b)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The solution for $\\bar{x}$ will be a combination of sine and cosine functions at\n",
"given natural frequencies, depending upon initial conditions, substituting\n",
"$\\bar{x}=\\bar{u}\\sin(\\omega t)$ Eqn 4b becomes\n",
"\n",
"$-\\omega^2\\left[\\begin{array}{cc}\n",
"m_1 & 0 \\\\\n",
"0 & m_2 \\end{array}\\right]\n",
"\\left[\\begin{array}{c}\n",
"u_1 \\\\\n",
"u_2 \\end{array}\\right] \\sin\\omega t=\n",
"-\\left[\\begin{array}{cc}\n",
"k_1+k_2 & -k_2 \\\\\n",
"-k_2 & k_2+k_3 \\end{array}\\right]\n",
"\\left[\\begin{array}{c}\n",
"u_1 \\\\\n",
"u_2 \\end{array}\\right]\\sin\\omega t$ (5)\n",
"\n",
"where $u_1$ and $u_2$ are amplitudes of the sine function and $\\omega$ is the\n",
"natural frequency. Now, the two ordinary differential equations have been\n",
"converted to one eigenvalue problem. The eigenvalues, of $\\bar{\\bar{M}}$ and\n",
"$\\bar{\\bar{K}}$ are equal to the natural frequencies squared. In the following\n",
"example, the natural frequencies for $m_1=m_2$ = 0.2 kg and $k_1=k_2=k_3$ = 500\n",
"N/m. "
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"1st natural frequency is 16 Hz, mode shape: 1*x1(t)=1*x2(t)\n",
"2nd natural frequency is 28 Hz, mode shape: -1*x1(t)=1*x2(t)\n"
]
}
],
"source": [
"m1=m2=0.1 # 0.2 kg\n",
"k1=k2=k3=2*500 # 500 N/m\n",
"\n",
"M=np.array([[m1,0],[0,m2]])\n",
"K=np.array([[k1+k2,-k2],[-k2,k2+k3]])\n",
"e,v=linalg.eig(K,M)\n",
"w1=np.sqrt(e[0].real)/2/pi\n",
"v1=v[:,0]/max(v[:,0])\n",
"\n",
"w2=np.sqrt(e[1].real)/2/pi\n",
"v2=v[:,1]/max(v[:,1])\n",
"print('1st natural frequency is %1.0f Hz, \\\n",
" mode shape: %1.0f*x1(t)=%1.0f*x2(t)'%(w1,v1[0],v1[1]))\n",
"print('2nd natural frequency is %1.0f Hz, \\\n",
" mode shape: %1.0f*x1(t)=%1.0f*x2(t)'%(w2,v2[0],v2[1]))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"![Figure 3: Vibration modes of 2-mass system when masses are 200 g and 3 springs have stiffness k=500 N/m](./eigenvalues.gif)\n",
"\n",
"*Figure 3: Vibration modes of 2-mass system when masses are 200 g and 3 springs have stiffness k=500 N/m. The first and second modes are animated on the left and right, respectively.*"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The amplitudes of the displacements will be proportional to the initial\n",
"velocities. A **modal analysis** is concerned with the natural frequencies and\n",
"mode shapes. This example has two modes (because it has two degrees of freedom).\n",
"In mode one, the displacement of mass 1 is equal to mass 2 and $\\omega_{1}=\\sqrt{(k_1+k_2)/(m_1+m_2)}$. In mode two, the\n",
"displacement of mass 1 is equal and opposite to mass 2 and $\\omega_{2}=\\sqrt{2(k_1+k_2+k_3)/(m_1+m_2)}$. The natural frequencies\n",
"are 8 and 14 Hz, for mode 1 and 2, respectively."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### FEA modal analysis\n",
"\n",
"In the FEA modal analysis, each element is composed of 8 nodes (for hexahedrals) that are shared\n",
"by adjoining elements. Each node can move in three dimensions, which means the\n",
"total degrees of freedom of the system is $3\\times$(number of nodes). So the\n",
"eigenvalue solution would reveal $3\\times$(number of nodes) natural frequencies\n",
"and vibration modes. We cannot use the high vibration mode results because they\n",
"will not be accurate, but the lowest vibration modes will become more accurate\n",
"as we increase the number of nodes in the system. This is called \"convergence\".\n",
"The ideal FEA solution will reproduce our analytical or experimental results.\n",
"Table 1 demonstrates the convergence of the first natural frequency for a beam\n",
"where E=200 GPa, L=300 mm, t=3 mm, w=12 mm, and $\\rho$=8050 kg/m^3. \n",
"\n",
"*Table 1: Convergence of the first natural frequency for a steel beam \n",
"where E=200 GPa, L=300 mm, t=3 mm, w=12 mm, and $\\rho$=8050 kg/m^3. Error is calculated based upon best estimate $f_{best}$ and analytical prediction $f_{an}$.* \n",
"\n",
"|# of nodes | first natural Frequency (Hz) | relative error $|f_{best}-f_{prev}|/f_{best}$ | absolute error $|f_{an}-f_{best}|/f_{an}$|\n",
"|--- | --- | --- |---|\n",
"|1153| 26.252 | | 0.020\n",
"|1438 | 26.243 | 3.4e-4 |0.021 |\n",
"|1913 | 26.235 | 3.1e-4 |0.021|\n",
"|3621 | 26.219 | 3.1e-4 |0.021|\n",
"|4521 | 26.218 | 6.5e-4 |0.0217|\n",
"|10,437|26.217 |3.8e-5 |0.0218|\n",
"|21,351 | 26.215 |7.6e-5|0.0218|"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Text(0,0.5,'relative error')"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAawAAAEeCAYAAADWyiHSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3X1wFPeZJ/DvM6N3kEYSLzZICtZLjGMHYwvFie3dq4uDjF27uGwskVzV1d5dEkP2breusjZgkqq87F2ZAHFytXubO3Aqd3W1dVlbmPhCbtdEMr5U1uTFkjDGsY2NBObdgKSRBHrXPPdH90g9w0ia6emZnp75fqqmJPX8puehR/SjX/fv9/xEVUFERJTpfG4HQEREFA8mLCIi8gQmLCIi8gQmLCIi8gQmLCIi8gQmLCIi8gQmLCIi8gQmLCIi8gQmLCIi8oQ8twPIBiKyEcDG0tLSp26//Xa3wyEi8pSurq5rqrpsoXbC0kzOaWpq0s7OTrfDICLyFBHpUtWmhdrxkiAREXkCExYREXkCExYREXkCExYREXkCExYREXkCh7UTzeOVYxew9/BJXAyOYmV5MbZtWI3H761yOyyinMSERTSHV45dwM6DJzA6OQ0AuBAcxc6DJwCASYvIBbwkSDSHvYdPziSrsNHJaew9fNKliIhyGxMWUQxjk9O4EByN+dzFObYTUWoxYTlARDaKyP7BwUG3QyEHvH0+iD/5m1/P+fzK8uI0RkNEYUxYDlDVQ6q6JRAIuB0KJWFyOoQftn+AJ350FD1Xb8RsU5Dnw7YNq9McGREBHHRBBAA4dWUYX3/xOE5cmO0llxT48ad3r8Ch4xcxOhkCALSuq+KACyKXMGFRTguFFP/j6BnsefV9jE+FZrZ/5rYKPN96Dz6xpAQViwqw71e9AIDlpbwcSOQWJizKWecHRvBM23H8trd/ZluB34enH74dX/3jOvh9AgCoriiJeA0RuYMJi3KOquJA13l899C7uD4+NbP9UyvK8MMvrsUdt5ZFtK+pmO1VnWPCInINExbllGvXx7Hz4Am0v/vxzDafAH/+L+vxH79wOwrybh6HFNnD4pB2IrcwYVHOePWdy/jmz06g78bEzLbapYvw/da1WLeqYs7XVVt6WJcGxzA1HUKenwNsidKNCYuy3tDYJL7z8z/gYPeFiO1/dv8qPPvoHSgpmP+/QVG+H8tKC3F1eBzTIcWlwTHUVJbM+xoich4TFmW1N05dw7a247g4ODaz7dayIuxpuRv/4vZlce+npqIYV4fHARiXBZmwiNKPCYuy0ujENHa/+j7+59EzEdsfv2clvvvYpxEoyU9of9UVJeg+GwRgDLy4H0ucCpWI4sSERVnnrXNB/NVLb6HXUq2ioiQf//nxNfiTu1fY2qf1PhYHXhC5gwmLssbkdAh/e+QU/u71U5gO6cz2L9yxHLueXIPlpUW29229BHi+n0PbidzAhOUAEdkIYGNDQ4PboeSsDz8extdfegvvXBia2baowI9vbbwTm5tqICJJ7Z89LCL3MWE5QFUPATjU1NT0lNux5JpQSPGTN05jz+GTmLCUVrrvtko8v3mtY4MjaixzsTh5mMgdTFjkWef6jdJKvzsdWVpp24bV+PIf1c6UVnLCivIiiACqwOWhMUxMhWJOMiai1GHCIs9RVbR1nsdf/yKytNJdK8vwg833YPWtpY6/Z2GeH7eWFeHS4BhUgUuDo1i1ZJHj70NEc2PCIk+5OjyOnQffRsd7V2a2+QT4D59vwF8+9MmU9nqqK4pxyZzPda6fCYso3ZiwyDP+6cQlfPOVd9AfVVrp+c1r0fiJuUsrOaW6ogRvnhkAwKrtRG5gwqKMNzg6ie/+/A84eCyytNK/uX8Vnn30Uygu8KclDlZtJ3IXExZltH/+8Bq2HTg+cykOAFYEirC3ZS3+6JNL0xoLq7YTuYsJizLSXKWVNt1bhW8/dhcCxYmVVnJCdaWlh8XJw0Rpx4RFGefY2QE8/dJx9F6LLK303BNr8Ogae6WVnFDDHhaRq5iwKGNMTIXwt0c+xN+9fgqWykpY/6nleG5TcqWVnLAiUAS/TzAdUlwZHsfY5DSK8tNz/4yImLAoQ3zw8TC+/uJb+MPF2dJKiwvz8K0/vROtTdVJl1ZyQp7fh1vLinAhaPSuLgRHUb9ssctREeUOJixy1XRI8ZN/Po29v4wsrfTZ2kp8v9W50kpOqaksnklY5weYsIjSiQmLXHOufwRPtx3H762llfJ82L5hNb78YC18DpZWcooxUtCIlwMviNKLCYvSTlXxUuc5/PWhd3FjYnpm+6eryvDDzffgk7c4X1rJKazaTuQeJixKqyvDY3j25RM48v5saSW/T8zSSg3I92d2QVlWbSdyDxOWA7geVnz+8cQlfPNnJzAwMjmzrW7ZIvxg8z24p6bcxcjixx4WkXuYsBzA9bDmNzgyiW///B288tbFiO3/9oHbsOORO9JWWskJ1kEgF9jDIkorJixKqV9/eBXb2t7G5aHZ0korA0XY27oWDzakt7SSE24pK0K+XzA5rbh2fQIjE1MoKeB/I6J04P80SomRiSl875/ex//6zUcR2zc1VuE7j92FsqL0l1Zygt8nWFlejI/6jN7VhYHRjB4kQpRNmLDIcd1maaXTltJKlYsK8NwTa/DIp291MTJnVFfMJqxzAyNMWERpwoRFjpmYCuFvXvsQP/p/kaWVmu+8Bc89sQbLSgvdC85B1eUlAPoAcOAFUToxYZEj3r88hL968TjevRRZWunbG+9Ey7rMKK3klBpWbSdyBRMWJWU6pPjxr3vx/C8/wMT0bGmlz9UZpZWsa0hlC66LReQOJiyXvXLsAvYePomLwVGsLC/Gtg2r8fi9VW6HFZezfSN4uu2tmWXjAaO00o5H7sC/e+C2jCyt5ARrD4sJiyh9mLBc9MqxC9h58ARGJ43yRBeCo9h58AQAZHTSUlX8w5vn8J9+8S5GLKWV1lQF8IPNa7N+EEI1q10QuYIJy0V7D5+cSVZho5PT+M6hP+AztZWoKi+e45XuuTI0hh0vv43XT16d2eb3Cf7i8w34Cw+UVnLCssWFKMjzYWIqhODIJIbHJlHq0WH6RF6SUMISkQCMUtXtqvpIakLKHReDsS8nBUcm8eD3jmDVkhI8UL8UD9QvwQP1S7Bksbuj7P7v25fwzVdOIGgprVRvllZa65HSSk7w+QTV5cUzKyKfHxjFp1YwYRGlWkIJS1UHReQMgPrUhJNbVpbPrq0Uy0d9I/io7yx++vuzAIA7bi3FA/VL8WDDEtxXW5m2v+oHRybxrZ+/g/8TVVrpyw/WYvsjq3Ny1d2qiuiEVeZyRETZz84lwVYAHSLytKo+73RAuWTbhtUR97AAIM8nqF+2CB/1j2BsMhTR/v3Lw3j/8jB+8sZp+H2Cu6sDeNDsgTWuqkhJ4vjVB1ex/cBxfDw0PrOtqrwYe1vvxgP13iut5BRrTUEObSdKDzsJqx9G0moTkS8B6ADw5lyNVfWgzdiyXnhgRaxRguNT0zh2NoijPX04euoa3joXxJRlNu50SHHsbBDHzgbxX18/hYI8H5pWVeDBBiOBrakKIC+J+0kjE1N47h/fw9//9mzE9pZ11fjWxjs9W1rJKazaTpR+oqoLt7K+QCQEQAFYxyzH2okAUFXNmetFTU1N2tnZmZJ9Xx+fwptn+nH01DUc7enDu5eGMN9HV1qYh8/WVeJ+8xLi6ltK45682/VRP55+6TjO9M32HJYsKsBzm9Zgw13eL63khEPHL+Ivf3oMAPDwnbdg/581uRwRkXeJSJeqLvifyE4P61nETlCUQosL8/D51cvx+dXLAQADNybw294+vNFzDUdP9c3cTwkbHp9Cx3tX0PGesVDikkUFuL9+yUwP7BOVJTclsImpEP5Lxwf477/qiSit9PCdt+C5TWuw1OVBH5nE2sM6xx4WUVoknLBUdU8qAqHEVCwqwKNrVuDRNSsAAJcGR3H0VJ9xCbHnGi4NjkW077sxgV+8fQm/ePsSAOM+1IMNS1CQ50PHu1fw8dAY/D6JuOxYWpiH7zx2FzY1VmVVaSUnWO9hnedcLKK0cGweloiUqerQwi0pFVYEivHkumo8ua4aqoozfSN449Q1/MZMYNZVfgFjkvJLnecjtlmT1QP1S7C3dW1GzgXLBEsWFaAo34exyRCGx6YwODKJQElu39cjSjXbCUtE7gGwG0ATgHJzWxDGAIzdqvq6IxFSwkQEtUsXoXbpIvzrz61CKKR4//IwjvZcwxunruH3p/txY2J6ztcHivPw91/5bNaWVnKCiKC6ogSnrlwHYFS8CJQEXI6KKLvZSlgisgvAdswOvOg1v9YBeBhAs4jsU9V/n3yIlCyfT3DnyjLcubIMX/3jOkxOh/D2+SCe/G+/idl+aHSKySoONRXFMwnr/MAIPl3FhEWUSgmPexaRJwHsADAIoFVVfaraYD58ADYDGAKwVUSecDbczCQiG0Vk/+DgoNuhxCXf78O6VXOXflrJy4BxYdV2ovSyM1FnK4xRgo2q+nL0k6p6AMA6GL2vrcmF5w2qekhVtwQC3voLe9uG1SiOmmxcnO/Htg2rXYrIW1i1nSi97FwSbALQrapn5mqgqr0i0g3gM3YDo9Sbb+IyLSyiajurXRClHKu157jH761igrKphpcEidLKziXBLgCNIrJqrgZmVfdGAKkp+0CUASInD48g0aoxRJQYOwlrH4z7U+0isjb6SXO4eyeM+1xtyYVHlLnKS/KxuNC4SDEyMX3TXDciclbCCcscVPEygAYA3SLSJyKHzUcfjB5YPYADqvpjZ8MlyhzGXCxLL4v3sYhSylY5b1VthTF8/QyACgDN5qMCwGkYw92/6FCMRBmLVduJ0sf2oAuzp3UAAESk1tx22qG4iDwhYqQgawoSpVRCPSwRCYjItIi8at2uqqeZrCgXRfawmLCIUimhhKWqgzAuA9anJBoij4ms2s5LgkSpZOceViuAJSLytNPBEHkNB10QpY+de1j9MJJWm4h8CUAHjArtManqQZuxEWW86HqCqsq1w4hSxE7C6oUxx0pg1Axch9grEIu53R/jOaKsECjOR1lRHobGpjA+FcLV6+NYXlrkdlhEWclOwnoWsRMUUU6qrijBu5eMtUvPD4wyYRGlSMIJS1X3pCIQIq+qqSyeSVjn+kfQ+IkKlyMiyk6ODGsnymVcF4soPTisnShJNax2QZQWHNZOlKTIHhaHthOlCoe1EyWJk4eJ0oPD2omSVGW5JHhhYBShkMLn41wsIqdxWDtRkhYX5qGiJB8DI5OYmA7hyvA4bg1waDuR0zisncgBNZUlGBgZBGBUbWfCInKerfWwrESkTERuE5EyJwIi8iJWbSdKPdsJS0SeMVcYHgDQA2CLuf2wiHzFofiIPKHGOlKwnwMviFLBVsISkcMAdmN2hWHrHeb7AOzn5GLKJRFV29nDIkqJhBOWiDwFoBnAMQB1qtoQ1aQWwBEAzexpUa6o5tB2opSz08PaCmOU4EOqeib6SVUNqmozgEEAX0suPCJvqGEPiyjl7CSsRgDdqjq0QLtOAHU29k/kOVXlsz2sS8ExTE2HXIyGKDvZSVi9cbbLmWQlIhtFZP/g4KDboZBLigv8WLq4EAAwFVJcHhpzOSKi7GMnYR0D0Cgiq+ZqICK1MBJWp93AvERVD6nqlkAg4HYo5KJqFsElSik7CWsXjFGB7SKy1rJdAUBE7oGRqBTGSEKinMCagkSplXDCUtVuGOWZGgB0m3OxFMA3zO+7YAx336OqR5wMliiTRQxt7+fACyKn2ZqHZZZnaoAxfF3MR4X59TUA61R1p1NBEnlBDRdyJEopO8VvAQCq2gtjPhYAYzVic4FHopzEycNEqZV0LcEwJivKddZ7WBfYwyJynGMJiyjXrSyfrdB+aXAUk5yLReQoJiwihxTm+XFLmTEXK6TGBGIicg4TFpGDIgde8D4WkZOYsIgcxIEXRKnDhEXkIE4eJkodrjhM5CBOHiZKHSdXHH7K3M4VhylncfIwUeo4ueJweNVhrjhMOavakrB4D4vIWVxxmMhBK8qL4DP/dPt4aBzjU9PuBkSURbjiMJGD8v0+rAjM3se6yLlYRI7hisNEDqviwAuilOCKw0QO48ALotTgisNEDuPkYaLU4IrDRA7j5GGi1OCKw0QO4+RhotTgisNEDrMmLPawiJzDFYeJHHZrWRHyfIKpkOLa9XGMTU6jKN/vdlhEnmdn4vBXY9UNZLIiMuT5fVhhWcyRy4wQOcPOJcH9AAZE5FURecLpgIiyQU1EiSZeFiRygp2E9QKMKhYPAzggIn0i8iNzdCARIeo+FgdeEDnCzijBrapaCeP+1Y9hDLT4GoAuEflQRJ4WkdscjZLIYzh5mMh5tpcXUdXXYiSvJQD2AugRkTdZ/JZyVXUlJw8TOc2RFYdjJK/TANYB2OfE/om8hj0sIufZHtYeTUQeAtAKYDOMOVkAsFCBXKKsVM2EReS4ZFYcLhORTSLyoohMA2iHsfSIwBhJ2Gz2uIhyzvLSQhT4jf9e/TcmcGN8yuWIiLwv4R6WiHwVRk9qfXgTgCCAlwC0qeprzoVH5E0+n6Cqohinr90AYPSyVt9a6nJURN5m55LgfvMrkxTRPKotCetc/wgTFlGS7CSsAwD2MUkRzS/yPhZHChIlK+GEpaqbUxEIUbaJXBeLAy+IkrVgwjJH/wFAp6oOmZOC415NmEuMUK6KXBeLPSyiZMXTw+qAsd5VM4zlRFphLMyocb4Hy1RTTuIyI0TOiidh7YWRnHrNn7sB7ElZRERZggs5EjlrwYSlqjuifn4NxiKNWUtEymFU63hKVQ+4HQ9507LFhSjM82F8KoShsSkMjk4iUJzvdlhEnmVnPaxdlvta87XbJCLP2AvLdW0Ayt0OgrxNRKIuC7KXRZQMO5UudgC4N452X4Jxr8tTRGQ7AFboIEdEDrzgfSyiZMQzSjCA2cnCYV8TkfvmeVk5jEoYwUSCEZEtAHpVtSOR1zlFRBphJNl6AD1uxEDZhfexiJwTz6CLShgjA8MUxgm9Po7XRie6OZn3jfaZr5kzYZlJZSdmS0N1ANilqt3xvtc82gBsVdVeEXFgd5TrWLWdyDnxDLo4LSLh5CQATsFILAuNFOxX1cEEYnlhoQYish5GkV1gdtRiC4AWEWlNZoCEiLQB6FbVuJMs0UJYtZ3IOXFVulDV0+HvRaQDQLt1m10i0gJjftdmLDDIweyBhZNVc/iyoSWJtYlIhaomdBnSEsd6ALWJvpZoPhx0QeQcO6WZHo6nnVnVfZ2q/vk8zXYCaIzzrbeYX/dY73GpaoeI7Def3wKz52feD1u3wD7bzV5ZM4yEORB1KbBNRKCqvD5ItkQPulBV8HIzkT22F3BcoERTOYBnYfRY5kxYqjqTUMxeTts8b9lsfn0xxnP7YCSrZpgJK8FLe7tx8+rIXTBGRLoyAISyQ0VJPkoK/BiZmMb18SkERyZRsajA7bCIPMnOelgBGCfxhXpGAmdP9nUAEGtwhap2m3+1xl3jMOr1vdHbzP31OjSYg3KUiKCmogQnPx4GYPSymLCI7LEzD2snjEttp2GUbTpmbt9rPg7CSFZt8V4+jFM8ychWwkqGiGwRkU4R6bx69Wq63548ILJqO+9jEdll55JgCwBV1QZgZpj5mwD+QVXfMrdtAbBLREpVddixaOef1xWEg9Up4r1vZV563A8ATU1N8RYEphzCqu1EzrDTw6qD5VKfeclMYBngYJ7EB+B8pYv5EhJLKVFGYtV2ImfYSVixejndAJpibGuO0daum+4z2WxDlFasdkHkDDsJqxc3J6dOzFaesHLynlIvMHMJMoJlW8JzsIhSjZOHiZxhJ2G9BKBCRF4VkXvMbd0A6kTkCQAQkVoY97qc7PGEh7x/McZz4W3RQ9OJXBddnkmVtzqJ7Eg4YanqHhgjAx+GMWIQMOZGCYADIvIhjPJNCgfvYVnmVW239rLM77dHtSHKGIGSfJQWGeObRien0XdjwuWIiLzJTg8rPOF3K4zeFsyagQ8DOAOjKO4ggL2q+mNnwpyx1fzaJSLtItIOY4Kv9TmijGO9LMj7WET22K50oaovRP3cgfgquNumqvtFpBdGzy18z6wbwA63liQBABHZCGBjQ0ODWyFQhqupKMZ7l4YAGJcF7/1EhcsREXmP7YTlNLOm34Jzn8zEtFCNwLRS1UMADjU1NT3ldiyUmTjwgih58Szg+FAyb6CqR5J5PVE2qKlktQuiZMXTw+qAMYDCLn8SryXKCuxhESUvnoS1F8klLKKcF1HtgoMuiGyJZ8XhHekIhCibRSSs4ChCIYXPx3WxiBJha1i7lYiUichtIlLmREBE2ai0KB/lJfkAgImpEK5eH3c5IiLvsZ2wROQZEemDUeS2B+aKwCJyWES+4lB8RFkjsuIFLwsSJcpWwhKRwzDmQlXAWBfLem3jPgD7ReTV5MPzBhHZKCL7BwcH3Q6FMhirthMlJ+GEJSJPwajCfgxAXXhdLItaAEcANOdKT0tVD6nqlkAg4HYolMGs62Kx2gVR4uz0sLbCGDX4kKqeiX5SVYOq2gyjPNPXkguPKHuwh0WUHDsJqxFAt6oOLdCuEy4sWU+UqSLWxeI9LKKE2V0PKx5MVkQW0cuMEFFi7CSsYwAaRWTVXA3M9bDqYPSyiAhAlaWHdTE4iukQ5+MTJcJOwtoFY1Rgu4istWxXADAXdeyEw+thEXldSUEeli4uAABMTis+HhpzOSIib7GzgGM3gGcBNADoNudiKYBvmN93wRjuvoeFb4kiVfGyIJFtdhdw3AMjYR2B0dsSGElKALwGYJ2q7px7D0S5qcY68IJD24kSkswCjr0w5mMBAEQkYK48TERzYNV2IvsS6mGJSEBEpmNVscjlZMVKFxQvDm0nsi+hhGUmpTMA6lMSjUex0gXFy1rtgvUEiRJj5x5WK4AlIvK008EQZbuIHlY/LwkSJcLOPax+GEmrTUS+BGNF4jfnaqyqB23GRpR1qspnE9bloTFMTYeQ5096lR+inGAnYfXCGMYuANaZj1gzIMXc7rcdHVGWKcr3Y3lpIa4Mj2M6pLg0OBZxmZCI5mYnYT2L2AmKiOJQU1mCK8PGAo7nBkaYsIjilHDCMudgEZFN1RXF6PpoAACHthMlghfPidIsYpkRTh4mihsTFlGasWo7kT1MWERpZq12wcnDRPFjwiJKs5pKrjxMZAcTlgNYmokSsSJQDBHj+8tDY5iYCrkbEJFHMGE5gKWZKBEFeT6sKCsCAKgaizkS0cKYsIhcwKrtRIljwiJyQXUlq7YTJYoJi8gFkT0sJiyieDBhEbmAVduJEseEReSCGvawiBLGhEXkgojyTBx0QRQXJiwiF6wIFMHvMyZjXRkex9jktMsREWU+JiwiF+T5fVgRKJr5+QLnYhEtiAmLyCXW+1jnWLWdaEFMWEQu4X0sosQwYRG5hFXbiRLDhOUAFr8lO1i1nSgxTFgOYPFbsiOi2gXvYREtiAmLyCXsYRElhgmLyCXLS4uQ7zfmYvXdmMDIxJTLERFlNiYsIpf4fYKqcvayiOLFhEXkIlZtJ4ofExaRi1i1nSh+TFhELqqpZA+LKF5MWEQuYg+LKH5MWEQuiriHFWQPi2g+TFhELqphPUGiuDFhEbloWWkhCvOM/4bBkUkMj026HBFR5mLCInKRiKCKvSyiuDBhEbmsmutiEcWFCYvIZbyPRRQfJiwil3FdLKL45LkdQDYQkY0ANjY0NLgdCnnQ5cHZXtVPf3cWa6vL8fi9VS5GlDteOXYBew+fxMXgKFaWF2PbhtUZe+ydiNVL/95YRFXdjiFrNDU1aWdnp9thkIe8cuwCdrz8NsanQjPb/D7B52orULtssYuRZb/TV6/jt6cHMB2aPQf6fYLP1lagdmlmHfvT167jd0nGGmsfxfl+7Nq0xvWkJSJdqtq0YDsmLOcwYVGiHvzeEVwI8r4VuaeqvBhvPPuQqzHEm7B4D4vIRReZrMhlXvod5D0sIhetLC+O2cMKFOfjmYdvdyGi3PH9X36AwdGbJ2oHivPxzIbVLkQ0t+8fPpl0rHPtY6VlTbZMx4RF5KJtG1Zj58ETGJ2cntlWnO/Hdx+7y/X7CtmutCjfM8e+tDAv6VhLC/Ow7cBxTE5H3sPalmHJeT5MWEQuCp9svDxyy6u8dOydiPXxe6vwRs81tHWeBwCUFPjx3BPuD7hIBAddOIiDLogok/3v353FN352AgDwr+6rwa5Nd7sckYGDLoiIKKswYRERkScwYRERkScwYRERkScwYRER5SAvjrfjKEEHichVAB+5HYeDAgAG3Q7ChkyIOx0xOP0eTu0vmf3YeW2ir1kK4FqC70Gp/Z1eparLFmrEhEVzEpH9qrrF7TgSlQlxpyMGp9/Dqf0lsx87r030NSLSGc8QaoqUCf+veEmQ5nPI7QBsyoS40xGD0+/h1P6S2Y+d12bC550LXD/O7GERUU5hD8u72MMiolyz3+0AyB4mLCLKKaoaV8ISkXIRGRCRllTHRPFhwiIiiq0NQLnbQdAsJiwioigish1ApdtxUCQmLCLyNBHZIiLrHdxfI4DdAFqd2ic5gwmLiDxLRMoB7MMCyUVEGkWkzbwnNWB+3zhH8zYAW1W11+l4KTlMWETkZS8s1MDsfXUBaAHQbz5aAHRFD6gQkTYA3fEOzKD04orDROQpZpJpBrAZCwyKMHtg7eaPzaraYW5fb25vE5EKVQ2a+10PoDZlwVNS2MMiIq/ZCWAL4hvBFy4ltCecrADA/H5/VJtmc58DIqIiEq6q0Gb5nlzEhEVEnqKq61RVVFWw8MCIZvPrizGe2xfVZjeAdVEPANhh+Z5cxEuCRJTN6gBAVbujn1DVbhGxtrlpkIX5fG+s11P6sYdFRNmszqE2lAHYwyKibBdc4Lk574WZlx0pQ7CHRUTZbr7BGSy95CFMWESUzeKZ/MsJwh7BhEVE2awXmCm3FMGybb5LhpRBmLCIKJu1mV+/GOO58LZ9MZ6jDMSERURZy1Jiabu1l2V+vz2qDWU4JiwiynZbza9dItIuIu0wagtanyMP4LB2Ispqqrp+OyaGAAAFlUlEQVRfRHphVLIIL0PSDWCHtVwTZT5RZYksIiLKfLwkSEREnsCERUREnsCERUREnsCERUREnsCERUREnsCERUREnsCERUREnsCEReQAEekRkayY1CgijSLSJSIqIhldZ8+MsWvhlpQNWOmCiKK1wViFtwOzJYyIXMeERUTR6gB0q2qz24EQWfGSIJFHiYjjq+Va9tnv9L6JksWERRlNRLab9ykazUe7+fOA+X1jVPt94fYx9rXefG73HPtvsdy76Ylqty98n8ps0zJPzLuj2u6ep22jiLRFtd8yz3GoM1/TA2AgjkMYfn2LebwGLMeuJarNPss+w8dq3ntYiX4+icQT1X635bPpEpHtC8QV13Gdo22P+Xk7/gcBJUlV+eAjYx8w1ixSy9ceGPdYesyfFUCdpf0+c1tjjH2tN5/bHWP/4de1m4/wvttg3McZsHyv0e9hiSf82oGotj0Ayuf4t6nZ1tq+fY62681990S3mecYts3zPrujjs9uS7zbAax38vNJJB6zbXnU8+HPwnqsu5I4ruuj2rZb9t/j9u8/H1G/b24HwAcf8z2iTj7bo54Ln/isJ127CSs6AUWfyMpjvId1P9YTdItlu/WEu8+yvc5ygq+bo/32GHEORB+HBY5fyxzvU2c5MTdGvf9NJ3YHP59E4wkf6+jPoMXyvl1R+0nkuIa3rY+KPZwMW+I5Dnyk5+F6AHzwMd/DckLsivFcOKm0WbbZ7mHFaH/TCXSe9+2J3rfluXAS0PBJ13JCjBVn+KQ7EM9xWOD49czzPuGTfnuMWBNNWPF+PnHHE+u4RbXfF/3eNo6rAtAYbcMrEt+0Hz7ce/AeFnnFizG29Tq4/1jDt/sBQFW7E3jfm+75qGoQxhBxwDhpAkATgGCMfUNVe833KI9xHyXRBQfr5nmfA5ZYkhXv55NIPOFjdcA8htFi3V9L9Lh2A4B5D229pW23qu6JtR9yDxMWeYWTySkWR0bFmSfFWMInvvBJuBzGiVNjPSztKqP282a8sYhIeB+d8zTrNWNJ1oKfj414mizb4n3PRI9rq7mf9QDCg0DaRWQLB11kHs7DIkozy4kwCGDXAs3TMrxcRMrn6MW4wjxG8/7bVTUoItGvARI4ruYfGPVm76oVRuIKP3aLyBfYy8ocTFiUS+oWbpIcEambo5cVHt7dbZ5ogwCgqntSFYuq9pon9Pku+YUv0aU8WSUaj4j0WrbdJHrIfDLHVVU7YF5uNXuCOwBsgTFwpD6RfVHq8JIgZavoS2mA8Rd0qm2N3mD+5b8exok4fBLuhHHpKtZ8sfLwpSkH4gnfs4n1PuF5T/NdonNaIvGEj1XLHJfnbjrWSOC4mnPaekSkzdpOVXtVNbzvlP+RQ/FjwqJs02N+jTiZmZNG19/c3HHbrRNgzRPta+aP1stUO8yvr1nu7YS1wbgX04bkhd+nzXrSN9/zBfPHOSc2p0Dc8Zi9vnBP6bWo9i0wekBz7X/B42r+8VAHIyFG/G5YPkNeDswgTFiUbcIjzVrCfz2LUc17H2ZPZqnSC+MEF66aEJ7k2gigw3qZyrwvsgfGCbTHbN8uIgMwEmuHqu5PNiBz5N0BmPOczIoPXTAnMgPYb14OSwsb8eyCcUwbLe3Dk5M7EDXwwsZxDf9OtFvahvcPAE8596+nZDFhUVYx/2peB+NkVgljbg8ANJsJYweMuTpO2w1jDtY6GCfMShgn2W4AOzRGIVlV3QGg2RJrE4wT8NZY7e1S1VYYPc4OGIkiXIm91XLpK20SiUdVg5ZjGk5clQD2mMdoN6KGtydyXM3fiVZz35WY7YXvB1DPAReZRdSYJEdERJTR2MMiIiJPYMIiIiJPYMIiIiJPYMIiIiJPYMIiIiJPYMIiIiJPYMIiIiJPYMIiIiJPYMIiIiJP+P8kzHsBe4fD9gAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"n=np.array([1153, 1438, 1913, 3621, 4521, 10437, 15316, 17801, 21351])\n",
"f=np.array([26.252,26.243,26.235,26.219, 26.218, 26.217,26.216, 26.215,26.215])\n",
"rel_error= np.abs(f[1:]-f[0:-1])/f[1:]\n",
"abs_error= np.abs(f-26.8)/26.8\n",
"plt.loglog(n[1:],rel_error,'o-')\n",
"plt.xlabel('number of nodes')\n",
"plt.ylabel('relative error')"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Text(0,0.5,'absolute error')"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEeCAYAAACAMA7hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8VOW9P/DPN3sCWQHZt4CCCgIJpdValwKite1VWWxd+uomYHu9tlVAWn/3XvtrS4PWbr7EoP3d9kprNUhpvV6lhKq11VaYhEVRloR9h0xCyEKSmef3xzkzOZlMljlzZs6ZM5/36zUvcs48c+bLSeabJ895zvcRpRSIiMh9UuwOgIiIYoMJnojIpZjgiYhcigmeiMilmOCJiFyKCZ6IyKWY4ImIXIoJnojIpZjgiYhcKs3ONx88eLAaN26cnSEQESUcj8dzVik1pK92tib4cePGYdu2bXaGQESUcETkUH/acYiGiMilmOCJiFyKCZ6IyKWY4ImIXIoJnojIpWydRUNElAg2Vh/D45v24Hh9C0YUZGPZvEm4bcZI24/VFyZ4IqJebKw+hpUbdqKl3Q8AOFbfguXrd2LvqUZcO3FwRMf62/6zeO7tA2jzdR5r5YZdABCTJC92Ltk3c+ZMxXnwROREPr/CtoN1+Oqvt6KpzRfT9xpZkI2/P/LpfrcXEY9SamZf7diDJyLS+fwK/zxwDq/tOonXPziJM40X4/K+x+tbYnJcJngiSmodPj/+eaAOr+46gT9/cBJnL7T163UZaSkoHVMY0Xt5DnvR1uHvtn9EQXZEx+kvJngicr3QC5vfmXsphuRm4bX3T2DTB6dQ1xQ+qQ8emIlJQwdg66H6Lok5Oz0Vq+6YGvG4uTaevwst7Z1DPtnpqVg2b5K5/1gfmOCJyNVCk+qx+hY8VLGzx/aX5GbilinDcMvU4fjYuCKkpohlM18Cr4nXLBpeZCUiV/D7FU43XsThumYcqWvGEW8zjtS14JUdx4OzVnoyNC8Tt0wZjs9MHY6ZYwuRkiJxitocXmQlIldRSqGhpR1H6lpwxNtsSOQtOFrXjKP1LWHHt3vz1U+Ox61XDcOM0c5P6mYwwRORY7S0+XDU29n7DvTED9dpSbzxYodl7zWiIAv//rkrLDueEzHBE1HcdPj8ONHQqifw5mBvPNATj3ZaYmFOOkYX5WiPwhyMLsrGkbpm/NffD+JiyEXS5fMmR/vfcTwmeCKKmvEi5NC8LHxh1miMHzwAR71aL/yw3hM/Ud+KDr/5637Z6akYXZStJ+8cjCrMxphAQi/KwcDM8Clt8rC8uF3YdBJeZCWiLlrbfWhs7UBjazvO6/8Gt1uM+7Wva85cQO3ZJliRSlJTBCMKsrSkbUjio4tyMKYoB4MGZEDEfWPlkeJFVqIk1O7zBxNvY2sHzrd0TdLnDcm663ZnAu9rxkm0huRmYrQhaY8uzMEovVc+PD8LaakscmsVJngih/D5FS4Ykm63ZNzSjsaLXXvQ51vaDT3sji430DjB168drw+faEMpowpzkJWeandYSYMJnsgCfr9CU1tH195wMFF3Hd4IP/zRgQsWzhCJRnqqIDcrHblZacjNSkNe8Ov0LtuBf//PH98Pe3v/yIJsPPpZd89ScTomeEoaPd2NqJRCS8i4s9Yz7p6sjUMZxh72hYsdiOLaoWVSBIbknI68Lv+mIS+7a7Lu3iYdWekpEY1zX+zwx/X2e+o/XmQl12v3+bHmzf345V/2o93X+fMuALLTU9DmU1HN7LBSbmZa195ydmdPWkvCXXvWedldk/WAjFRbLkLGcxEL4kVWSmJ1TW2oOuSF57AXVYe82HG0Hq3t3S8cKgDNYfablZ2e2i3h9tiDzgxN4OkYmJmG1AS9m/K2GSOZ0B2ICZ4Sms+vsPdUI6oOe+E55EX14XocONsU8XEy0lK0HrGxt5zdmYjD9ajzDD3qgVlpSOfsD3IYJnhKKA0t7ag+7EXV4XpUHfJi+5H6fl2cTBWBL8xw5LD8LLy17AZkpnFmB7kPEzw5Stex3Cx86epxKByQgapDXlQd9mLf6Qt93lCTkZqCKSPzUDKmEKVjC1EythDv1pwLeyHwkZsnM7mTazHBk2NsrD6GRzbsDI6XH6tvxarXPurzdZfkZmqJfIyWzKeMzOuWtONdh5vICZjgyTFWb/oo7MVQo9QUwRXD81A6thAzxhSgdGwhRhZk92vmCC8EUrJhgifHOF7f2uNzy2+ehJIxhbhqVD5yMvhjS9Qf/KSQI/xpx/EenxtZkI1v3DAxjtEQuQPndZHt9pxsxIr14dfI5B2RROaxB0+2amhpx5LntwVntwwZmIG01BScbGjlhVCiKDHBk238foXvvLgdB881AwByMlLx2/s+gcuG5tocGZE7cIiGbPPUG/ux5aPTwe3HF0xjcieyEBM82eKNPafx08q9we3F1xXj1quG2xgRkfswwVPcHTrXhAdfqA7ekXp18SAs54VUIssxwVNctbT5sOR5D863avVjhudn4Zd3zeAybUQxwE8VxY1SCis37MRHJxsBaDVj1txTisEDM22OjMidmOApbn7zzkFs3N55Q9P3/+VKTB9dYGNERO7GBE9x8d6BOvzg1Q+D21/42Gh8YdYYGyMicr+IEryI5IuIT0Rej1VA5D6nzrfiG7+tCi6LN21UPv7z81faHBWR+0WU4JVSDQAOApgQk2jIddo6/PjGb6tw9sJFAEDRgAysuacUWemswU4Ua2aGaBYCGCQiD1kdDLnPD17dDc8hLwAgRYCnvjgDIwqybY6KKDmYKVVQBy3JV4jIFwBUAtjaU2Ol1AaTsVGCe9lzFP/97qHg9iO3TMY1EwfbGBFRcjGT4GuhLUgvAEr1R7hF1ETfz7/Fk9D7xxrw3T/sCm7fOnU47vtUsY0RESUfMwn+EYRP6EQAAG9TG5au8+Bih7Y606WXDETZgqv6teoSEVkn4gSvlFodi0DIHXx+hQdf3I6j3hYAwMDMNDxzbykGZrJwKVG8WTYPXkTyrDoWJa6fbt6Lv+49E9z+yaJpmDBkoI0RESUv0wleRKaLyCYROSciPgBe/evXReRGC2OkBPHnD07iqTf2B7f/9caJmHflMBsjIkpuphK8iKwC4AEwF0AhgAP6oxDATQAqReRpq4Ik56s5cwHfeWlHcPu6y4bg23MvszEiIoo4wYvIfAArADQAWKiUSlFKTdQfKQAWATgPYImI3G5tuORETRc7sPR5Dy5c1CpEjirMxs/vnI7UFF5UJbKTmR78EmizaEqUUi+HPqmUWg9t6qTobcnFlFJYvn4n9p2+AADITEvBM/eUonBAhs2REZGZBD8TQJVS6mBPDZRStQCqAHzMZFyUIJ59uxav7joR3P7R7VMxZWS+jRERUQCrSZJp7+w/ix+/9lFw+0tXj8X80lE2RkRERmYSvAdAiYiM7amBiOQDKAGwzWxg5GzH61vwry9UQy8QidKxhXj01ivsDYqIujCT4Muhja9vFpFpoU+KyHRoiV0BqIguPHKi1nYf7l/nQV1TGwBg8MBMPH13CTLS+AchkZOYuZN1vYi8DGA+gCoRqUdnT30mgAJovwAqlFLPWRYpOcZjr3yAHUcbAABpKYKn7y7B0Lwsm6MiolCmulxKqYXQpkMehDb3fS66zolfqJS606IYyUFeeO8wXnjvSHD7e7dejlnji2yMiIh6YrpAiD4dcj0AiMh4fd8Bi+IiB9p+pB7/8ccPgtu3TR+BL18zzr6AiKhXlizZp5Q6wOTubmcvXMT96zxo82kVIicPy8WqO1ghksjJuGQf9anD58cDv6vGiYZWAEBeVhrK7y1FdgZL/RM5GZfsoz49vmkP3q09BwAQAX7+xRkYO2iAzVERUV+4ZB/16tWdJ1D+19rg9rdmX4YbJ11iY0RE1F9cso96tPdUI5at76wQOXvyJXjg0xNtjIiIIsEl+yis863tWPq8B81tPgDAuEE5ePLO6UhhhUiihMEl+6gbv1/hoZd2oPZsEwAgOz0V5ffORH52us2REVEkLJkmSe6y5q0abN59KrhdtuAqTBqWa2NERGQGp0lSF2/tPYMn/rwnuP21a8fj89NG2BgREZnFaZIUdKSuGf/2QjWUfoXl4+OL8Mgtk+0NiohM4zRJAgC0tPmw5HkPGlraAQDD8rLw1F0lSE9lhUiiRMVpkgSlFL63cRd2nzgPAEhPFTx9TwmG5GbaHBkRRYPTJAnr/nEIG6qOBbf/8/NXomRMoY0REZEVOE0yyXkO1eGxV3YHtxeWjsJds8bYGBERWSXqAVYRyRORcSKSZ0VAFD+nG1tx/7oqdOjr7k0dmY//e9sUVogkcgnTCV5EHhaRcwC8AGoALNb3bxKRr1kUH8VIu8+Pb/62CqcbLwIACnPSseaeEmSl85IJkVuYSvAisglAGTpXcDJ2+WYBWMuboZzth69+iK0HvQCAFAF++cUSjCrMsTkqIrJSxAleRO6DtjxfNYBipVRo9anxAP4CYC578s70h+qj+PU7B4Pby+ZNxrWXDrYvICKKCTM9+CXQZtF8Wil1MPRJpVS9UmougAYAS6MLj6y2+/h5rNywK7h9y5RhWHp9sY0REVGsmEnwJQCqlFLn+2i3DQAzh4M0NLdj6ToPWtu1ZfcmDBmAxxdO40VVIpcyk+Br+24CgMndUfx+hQdfrMbhumYAwMDMNJTfOxMDM02vu05EDmcmwVcDKBGRsT01EJHx0BL8NrOBkbV+tmUf3txzJrj9xMKrMPGSgTZGRESxZibBr4I2a2aziEwz7FcAICLToSV2BW2mDdmscvcp/GLLvuD2/TdMwM1ThtsYERHFQ8QJXilVBa1cwUQAVfpceAXgu/rXHmjTJ1crpf5iZbAUuQNnm/Dtl7YHt6+dOBgP3zTJxoiIKF5MzYPXyxVMhDYdUvRHof7vFgClSqmVVgVJ5jS3dWDp8x40tnYAAEYWZOMXX5yBVC67R5QUTF9hU0rVQpsPD0Bb7UlfEIQcQCmFFS/vwp5TjQCAjLQUPHNPKYoGZNgcGRHFi2XFvpncneVXfzuAV3YcD27/4LYpmDoq38aIiCjeuJqDC71bcw6rXvsouH33x8dg0czRNkZERHZggneZEw0teOCFKvj0CpEzxhTg3z93hc1REZEdmOBd5GKHD/evq8LZC20AgMEDM7Dm7lJkprFCJFEyYoJ3kcde2Y3tR+oBAKkpgqfuKsGw/CyboyIiuzDBu8RLW4/gd/88HNz+7mcuxyeKB9kYERHZjQneBXYercejf3w/uP25aSPw1U+Osy8gInIELtmX4Oqa2nD/uiq0dWgVIicNzUXZ/KmsEElEli7Zd5++n0v2xUmHz48HXqjCsfoWAEBuVhrK7y1FTgYrRBKRtUv2BbqMXLIvTp748178ff+54PbP7pyOcYMH2BgRETkJl+xLUK/tOoFn3qoJbj84+1LMvnyojRERkdNwyb4EtP90Ix6u2BHcvnHSEDw4+1IbIyIiJ+KSfQmmsbUdS573oKnNBwAYU5SDn905AymsEElEIbhkXwJRSuHhih2oOdMEAMhKT0H5vaXIz0m3OTIiciIu2ZdA1rxVg00fnApu//iOq3D5cM5OJaLwuGRfgnh73xk8sWlPcPsrnxyH22aMtDEiInI6LtmXAI56m/FvL1RDLxCJWeOK8N3PXG5vUETkeFyyz+Fa231Yus4Db3M7AOCS3Ew8dfcMpKeyygQR9Y5L9jmYUgqPbnwf7x/TJiylpwrW3FOKS3JZIZKI+mbmRqdVIvLp0P2hyV1E5ovIw9EEl+x++8/DWO85Gtz+989egdKxhTZGRESJxMzf+SsAzOhHuzvBi6ymVR324rFXPghuzy8ZhXs+0ePEJSKibvocohGRfABrQ3YvFZFZvbysAMAcAPVRxJa0zjRexP3rPGj3aVdVrxyRhx/ePoUVIokoIv0Zgy8CsNCwrQBM0B99Cf3FQH1o9/nxzd9V4dT5iwCAgpx0PHNPKbLSueweEUWmzwSvlDogIoFkLgD2AygHsLqPl9bxomvkfvzaR3jvQB0AQAT4xRdmYHRRjs1REVEi6tcsGqXUgcDXIlIJYLNxH0VnY/UxPL5pT7Cue8DDN03CdZcNsSkqIkp0EU+TVErdFItAktXG6mNYuWEXWtp9XfZPGZGH+6/vzygYEVF4ESd4EVkVQfPNvJu1d49v2tMtuQPAuaY2VogkoqiYudFpRR/PK8PX56Dd7Uo9OB4yLBNwsqE1zpEQkduYSfBze3muGMAiALMBLFdKPWEqqiRSOCADdU1t3faPKMi2IRoichMzY/Bbenl6C4BnRWQxgDUi4lFKvWE6Opc7e+EiLoYZnslOT8WyeZNsiIiI3CQmFauUUmuh1Y3nnaw9UEph5YZdwZWZUkSbgzqyIBur7pjKUsBEFDXTxcb6oRbA/BgeP6FVeI5i8+7OxTt+/ZVZnBJJRJaKZc3ZErBUQVhH6prx/Vd2B7e/dPVYJncispyZaZLT+2hSBGAJgPEAKs0E5WZ+v7au6oWLHQCA4sEDsPIWLt5BRNYzM0RTha5TIcMJTODua0pl0vl/fz+Af+qlCFJTBD9ZNA3ZGawzQ0TWM5Pgn0XfCb4GwHqWM+hq76lGrDasq/rNGyZgxhjWdyei2DAzTXJJLAJxu7YOP771++1o6/ADAKaMzMMDsy+1OSoicjMu7Bknv9iyD7tPaEvvZaSl4KeLpnNdVSKKqf4s+NFteb5IsBYN4DnkxdNv7g9ur7h5Mi4dmmtjRESUDPozRFOJvsfce5PUVxCb2zrw0Evb4dfP4NXFg/CVa8bZGhMRJYf+JPjHEV2CT2o/+t8PcfBcMwAgNzMNTyyaxiqRRBQX/VnRiVMdTXpr7xms+8fh4PZ/fP5KjGQRMSKKE8uu8olInlXHcoP65jYsq9gR3J535VDML2F9GSKKH9MJXkSmi8gmETknIj4AXv3r10XkRgtjTEiPbnwfpxu1hbMHD8zAj26fChEOzRBR/JhK8PqqTh5oteELARzQH4UAbgJQKSJPWxVkovnTjuP4n50ngts/vuMqDBqYaWNERJSMIk7wIjIfWgmCBgALlVIpSqmJ+iMF2oIf5wEsEZHbrQ3X+U42tOLRP+wKbt85czTmXDHUxoiIKFmZ6cEvgTarpkQp9XLok0qp9QBKodWjSaq7XpVSWLZ+B863aoXERhVm49HPspAYEdnDTIKfCaBKKXWwpwZKqVpoRck+ZjKuhLTuH4fw9r6zAAAR4MlF05GblW5zVESUrHivvEVqz1zAD//3w+D24k8VY9b4IhsjIqJkZybBewCUiMjYnhqISD60BT+2mQ0skXT4/PjOSzvQ2q4VEps0NBffnnuZzVERUbIzk+DLoY2vbxaRaaFP6guCbIM2Tl8RXXiJYc2bNdh+RFu8Kj1V8OSd05CVntQVGojIAcyUC14vIi9DW2+1SkTq0dlTnwmgANovgAql1HOWRepQ7x9rwM+37Atuf2vOZbhyRL6NERERaUyNwSulFkKbDnkQ2tz3ueg6J36hUupOi2J0rNZ2H7794nZ06JXESscWYun1E2yOiohIY2ZFJwDB6ZDrAUBExuv7kmoFp8c37cG+0xcAADkZqXhy0TSkspAYETmE6QRvlGyJHQDeqTmLX/2t87/9vVsvx9hBA2yMiIioq2hq0dyhX1ANbE8XkX0i4hOR98JdgHWL863tWFaxM7h946QhuGvWGBsjIiLqzkypgnwR2Q9thsycwD4AfwEwAdoF1pkAPCIyzrJIHeT7r+zGsfoWAEBBTjrK5l/FQmJE5DhmevBlAIoBbIG22hOgXXAtAFBuqEeTAmC5FUE6yaYPTmK952hw+4e3TcUleVk2RkREFJ6ZMfg5ALxKqZsM+xZCm/deBgSnUtZCm1njGmcaL2Llhs5CYrdNH4FbrxpuY0RERD0z04MvRvc7VGcCqA2pT1Ort3UFpRRWbtiJuqY2AMCwvCw89vkpNkdFRNQzMwm+FkCwyIqIzIY2PFMZ0q4YQL350JylYttRVH54Orj9xMJpyM9hITEici4zCb4aWi2awKpNKxBSlkBP+sXQfhkkvCN1zXjslQ+C21++ZhyuvXSwjREREfXNTIJfBW2mTKW+VN9saMMzfwEAEXkRwJ+hJf1VVgVqF59f4aGKHWhq8wEAiocMwIqbJ9scFRFR3yJO8EqpKmjL8m2Hluir9e2Amfr+FUqpDVYEaadf/a0W7x2oAwCkpgieXDQd2RksJEZEzmfqTlalVCW0VZvCmeOWO1v3nGzEE5v2Bre/eeNETB9dYGNERET9F1WpAhHJgzbWHpgtU+WW5N7W4ce3XtyONp9W4/2qUfl44NMTbY6KiKj/TCV4PbGXAVgc5jkvgFVKqZ9EGZutfr5lLz48cR4AkJmWgicXTUN6KhfAIqLEYbYHXwVgPPSLrfo2oA3bzAawWkTmKqVujj7E+PMcqsOaN2uC2ytunoyJl+TaGBERUeQiTvAi8mNoQzIeaOPtDSHPF0ArYzBXRB5KpJ78xupjKHv9I5xoaA3uu2bCIHz5mnH2BUVEZJKZMYc50KZAzg5N7gCglKqH1osXAAmz6MfG6mNYuWFXl+QOAHOvGIoU1ngnogRkJsGXQLuYer6nBnqSr0LPM20c5/FNe9DS7uu2/7m3XXHNmIiSkNlSBf2pMVOABCpVcFwv/9vf/URETmcmwb8MoFBEbu+pgYjMh/ZLYK3ZwOJtaA8lf0cUZMc5EiIia5i5k3UFtJkz60XkR8ZFPURknIisAvASAI9SaqVVgcba5cMHdtuXnZ6KZfMm2RANEVH0+pxFo9ebCfsUtEJjK8KsZiQACkTk60qp56ILMfYamtvx3gFvl30jC7KxbN4k3DZjpE1RERFFpz/TJA9CmzUTqYSZevL8Pw4Gi4ldNnQgXn/wOs6cIaKE12eCV0pNiEcgdmlt9+G//n4wuL30+glM7kTkCjG7915Evi4iT8fq+Fap2HYE5/RVmkYWZONz00bYHBERkTVMFxvTL672NF2yAMAj0MoZfMPse8Rah8+P8r92rkny9U+NZ70ZInINM6UK8qHNoinpqym6L+PnKK/uOoGjXm2ee2FOOu782GibIyIiso6ZHvxKaHeo1gJYD610wQwAj+vPTwAwH0CFUsqxpQqUUl0Kin35mvHIyYiqejIRkaOYyWgLACil1EQAEJESAFsB/F4ptV3ftxjAKhHJVUo1Whathd7cewYfndRCy05PxZeuHmtzRERE1jIz4FwMw9CLvoSfwFB3Rim1FoAXWs14RzL23r84awwKB2TYGA0RkfXMJPhw9WWqoK3FGrpvronjx5znkDe4zmpaiuDrnxpvc0RERNYzW2wsNJlvgzYWH6o/Rcni7pm3Onvv/zJ9JOvNEJErmUnwL0ErNva6iEzX91UBKA4UIBOR8dDG6mt7OIZt9p1qxObdp4LbS6935O8gIqKomSk2thpANYCboM2oAYAXoY3DrxeRfQD2Qytv4LgxeOO897lXDMWlQ7kUHxG5k6l5gUqpUhG5D0Cdvt0gIjcBKIc2TbIewFqnFRo7Xt+CjdXHgttLr3d1FQYiSnKmJ34rpZ4N2a6Eltwd61d/O4AOv1Y3bdb4IpSOLbQ5IiKi2Ema+/K9TW144b3Dwe372XsnIpdLmgT/3+8eQrNeEnjysFzcMGmIzREREcVWUiT45rYO/PqdzsWz779hAsIsUkJE5CpJkeBf2noE3uZ2AMCowmzcOnW4zREREcWe6xN8u8+PZ9/u7L0vvq4YaSwJTERJwPWZ7n92Hsexeq0k8KABGVhYypLARJQcXJ3g/f7QksDjkJ2RamNERETx4+oE/8ae09h76gIAYEBGKr509Th7AyIiiiNXJ3hjUbG7Pj4G+TnpNkZDRBRfrk3w2w7WYetBLwAgPVXwtWtZVIyIkotrE7yx9377jJEYlp9lYzRERPHnygS/52QjKj88DQAQARZfx7IERJR8XJngyw2995uuGIqJlwy0MRoiInu4LsEf9TbjTzuOB7dZEpiIkpXrEvxzb3eWBL66eBBmjGFJYCJKTq5K8HVNbfj91s6SwEtvYO+diJKXqxL8b945iNZ2PwDgiuF5uO7SwTZHRERkH9ck+KaLHfjNuweD2ywJTETJzjUJ/vdbj6BeLwk8pigHt0wZZnNERET2ckWCb+vw41dv1wa3WRKYiMglCf5PO47jeEMrAGDwwAwsKB1lc0RERPZL+ATv96suNzZ95ZPjkZXOksBERAmf4Ld8dBr7TmslgQdmpuGeT4y1OSIiImdI6ASvlMKaN/cHt+/++BjkZ7MkMBERkOAJfutBL6oO1wMAMlJT8NVrx9scERGRcyR0gjf23ueXjsTQPJYEJiIKSNgE/+GJ83hjzxkAWkng+z7FBT2IiIwSNsEbZ87cMmUYioewJDARkVGa3QFEamP1Max67UOcOn8xuI8lgYmIukuoHvzG6mNYuWFXl+SeIkDtmSYboyIicqaESvCPb9qDlnZfl31+pe0nIqKuEirBH69viWg/EVEyS6gEP6IgO6L9RETJLKES/LJ5k5AdUmcmOz0Vy+ZNsikiIiLnSqhZNLfNGAlAG3M/Xt+CEQXZWDZvUnA/ERF1SqgED2hJngmdiKhvCTVEQ0RE/ccET0TkUkzwREQuxQRPRORSTPBERC4lSin73lzkDIAmAGdtC8L58gE02B1EP9gVZyzf18pjR3sss68387pIXjMY/Pz2JlY/n2OVUkP6amRrggcAEdmmlJppaxAOJiJrlVKL7Y6jL3bFGcv3tfLY0R7L7OvNvC6S1/Dz2zu7P78conG+V+wOoJ/sijOW72vlsaM9ltnXm3ldovzMJQJbzyV78ERkGj+/zuaEHvxauwMgItP4+XUw2xO8UirsD4iIFIhIhYh4RUSJyGYRKYl3fETUs54+v0b6Z9krIgviERN1sj3B92ILgBIA9wEoDewTkQL7QiIiEyoA8HNrA0cWGxORYmjJvVQpVaXvnisiCsAcAOttC46I+k1ElgMosjuOZGVpD15EFovIHAsOVQdghSG5w9Bzr7fg+EQUwsLPb+B4JQDKACy06pgUGcsSvJ6Ay9HHN1NESgxj61796y5j60ogYJ+eAAAHaklEQVSpeqXUakP75dCGbCqVUpVWxUxEGis/vwYVAJYopWqtjpf6x8oe/LN9NdB7Bx4AC6D10uv0rz29XICZA2AJtCEb/qAQxYaln18RqQBQ1Z+LsBQ7USV4EVkgIuUi4oX2je6tbQGAzfrmXKXUBKXUBABz9X0V4S6gKqVW6+0mAFgkImXRxExEmlh9fvVkPwfaBAmyUbQ9+JUAFqN/V8gDt+uuNg6z6F+vNbbR/wzscnuv/mdeJbSePBFFLyafX2hJvwBAYIpz4G7KCsPXFAdRJXilVKlSSpRSgr4vpAR+078Y5rnykDZFAMrD9OiLwWEaIkvE8PNbBm1qs/EBACsMX1McxHOaZDEAGGfGBCilqkTE2KZSROqh/cYvgzbWFxiH5599RPEXyee3WydMf7423OspduJ5o1NxhG1KoU2JrIA2g6YYXefFE1H8RPr5JQeI941Ovc1hr4dhLFDvBXD+LJFz9PvzG0ofBqI4i3epgt4u5vBWZiJn4+c3wcQzwffn4igvoBI5Ez+/CSjuCT7cXW+GfSxDQORM/PwmoHgm+Ar93zvDPBfYVx7mOSKyHz+/CShuCd5wy/JyYy9A/3p5SBsichB+fhNTvC+yLtH/9egLeGyGVtvC+BwRORM/vwkmrtMklVJrRaQW2p1ugbKkVdBKA7NKJJGD8fObeGxfdJuIiGLDyUv2ERFRFJjgiYhcigmeiMilmOCJiFyKCZ6IyKWY4ImIXIoJnojIpZjgyRYiUuOW9Tn1NYQ9+vqjjq7Hosfo6bsluUG8F/wgcqMKaKsZVaLz1n0i2zHBE0WvGECVUmpuny2J4ohDNJQ0RMTyVYcMx6yz+thE0WKCdxkRWa6Ps5boj836tlf/uiSkfXmgfZhjzdGfK+vh+AsMY881Ie3KA+PsepsFvcRcFtK2rJe2JSJSEdJ+cS/noVh/TQ0Abz9OYeD1C/Tz5TWcuwUhbcoNxwycq17H4CP9/kQST0j7MsP3xiMiy/uIq1/ntYe2Nfr3m8v2OY1Sig8XPaDV5laGf2ugjRHX6NsKQLGhfbm+ryTMseboz5WFOX7gdZv1R+DYFdDGob2Gr1XoexjiCbzWG9K2BkBBD/83pbc1tt/cQ9s5+rFrQtv0cg4renmfspDzU2aIdzmAOVZ+fyKJR29bEPJ84HthPNeeKM7rnJC2mw3Hr7H755+PkJ83uwPgw+JvaNcP6/KQ5wKJwpikzCb40IQd+sEvCPMexuMYE9oCw35jgio37C82JMTiHtovDxOnN/Q89HH+FvTwPsWGRFYS8v7dEqGF359I4wmc69DvwQLD+3pCjhPJeQ3smxMSe+CXx4L+nAc+4vOwPQA+LP6GdiYQT5jnAkm4wrDPdA8+TPtuCaeX960JPbbhuUDSVIEkZUgg4eIMJClvf85DH+evppf3CSTJzWFijTTB9/f70+94wp23kPbloe9t4rwqACpM28DKTt2Ow4d9D47Bu9eLYfZZuep9uOmAdQCglKqK4H27jVkrpeqhTTkEtCQDADMB1Ic5NpRStfp7FIQZB450IYriXt5nvSGWaPX3+xNJPIFztV4/h6HCXR+I9LxWAYB+DWCOoW2VUmp1uOOQfZjg3cvKZB6OJbNG9CQSTiBRBJJWAbREo8I9DO2KQo6ztb+xiEjgGNt6aVarxxKtPr8/JuKZadjX3/eM9Lwu1I8zB0Dgou9mEVnMi6zOw3nw5HiGxFEPYFUfzeMyXVFECnroJdtCP0e9/t+VUvUiEvoaIILzqv9CnqD33hdCS/SBR5mIzGYv3jmY4Kk3xX03iY6IFPfQiw9MF6zSE1M9ACilVscqFqVUrZ4AexuCCQyZxDy5RxqPaOulBvZ1EzoFM5rzqrQ1WCv14xYDWAFgMbQLxRMiORbFDodoKCB0aAPQemixtiR0h96znAMtcQWS1jZoQwnh5usXBIYKLIgnMOYc7n0C8857GzKxWiTxBM7Vgh6GS7qda0RwXvV7CmpEpMLYTilVq5QKHDvmnQLqPyZ4qtH/7fLh129ymdO9ueWWG2/Y0RPTFn3TOGywQv93i2FsOqAC2lhyBaIXeJ8KY5LU3/NZfbPHG7FioN/x6H9VBHriW0LaL4DWw+7p+H2eV/2XbTG0XyBdfjYM30MOzzgIEzwFZmIsCPTORKs2WI7OD3+s1EJLCIG7IgM35ZQAqDQOG+jjuquhJZwavf1mEfFC+0VUqZRaG21A+syU9dDnmet3dHqg33gFYK0+PBEXJuJZBe2clhjaB26mqkTIhVYT5zXwM7HZ0DZwfAC4z7r/PUWLCT7J6b2yUmgf/iJoc6sBYK6eYFdAmytttTJoc+BLoSWYImhJqQrAChWmcJdSagWAuYZYZ0JLWEvCtTdLKbUQ2l80ldASa6BS5ELDUETcRBKPUqrecE4Dib4IwGr9HJUhZLpkJOdV/5lYqB+7CJ1/5a0FMIEXWJ1FlHaTAhERuQx78ERELsUET0TkUkzwREQuxQRPRORSTPBERC7FBE9E5FJM8ERELsUET0TkUkzwREQu9f8BmSFuXZ9k0VEAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plt.loglog(n,abs_error,'o-')\n",
"plt.xlabel('number of nodes')\n",
"plt.ylabel('absolute error')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"After 4,000 degrees of freedom, the first natural frequency mode has converged\n",
"(the relative error is stable and the absolute error has diminishing returns). The absolute error compares the analytical model\n",
"to the FEA model. There is a stable 2.2 % error in the frequencies predicted\n",
"from Euler-Lagrange beam theory and 3D FEA. **Make sure each reported frequency\n",
"has converged, demonstrate convergence in your report**. \n",
"\n",
"Check the accuracy of your model in 3 ways: \n",
"\n",
"1. Check convergence of your results (convergence)\n",
"\n",
"2. Compare FEA results to analytical model (verification)\n",
" \n",
"3. Compare FEA results to experimental results (validation)\n",
"\n",
"Method 1 proves that your results are repeatable. Method 2 is a verification of\n",
"your model. Verification of a model proves that with different methods, the\n",
"assumptions in your model will converge to the same solution. Method 3 is a\n",
"validation of your model. Validation uses measured results to confirm your model\n",
"is correct. For many engineering structures, method 2 is not feasible. Instead,\n",
"we check convergence with method 1 and validate the model with experimental\n",
"results. Here, we are using natural frequencies as our comparison metric. \n",
"\n",
"The other method we will use to confirm our results is to look at the mode shape\n",
"and the strain or acceleration at the point of measurement. We can probe the\n",
"value of a given measurement in the FEA results. The absolute value of\n",
"the result will depend upon the amplitude of the impulse applied (from the\n",
"hammer), but we can compare relative amplitudes based upon the frequency mode.\n",
"In the example used previously, 300-by-12-by-3 mm^3 steel beam (10,437 nodes),\n",
"the amplitude of axial strain on the top surface 11 mm from the support is given\n",
"for the first 3 vibration modes in Table 2. The mode 2 axial strain is nominally\n",
"zero. This is the first torsion vibration mode. In order to measure this response\n",
"we would need a different experimental setup. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"*Table 2: Relative amplitude of axial strain for first three modes on top\n",
"surface for a steel beam where E=200 GPa, L=300 mm, t=3 mm, w=12 mm, and\n",
"$\\rho$=8050 kg/m^3.* \n",
"\n",
"|mode | amplitude of strain (%)|\n",
"|--- | ---|\n",
"|1 | 1.2 |\n",
"|2 | 0.005 |\n",
"|3 | -6.5 |"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# References\n",
"\n",
"1. Fish, J. and Belytschko, T. A First Course in Finite Elements. Wiley 2007.\n",
"\n",
"2. [Ansys help topic: Modal Analysis](https://ansyshelp.ansys.com/account/secured?returnurl=/Views/Secured/corp/main_page.html)\n",
"\n",
"3. Engineering Vibration, D. J. Inman, Ch. 4 Multiple-Degree-of-Freedom Systems. Prentice Hall."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.5"
}
},
"nbformat": 4,
"nbformat_minor": 2
}