Axisymmetric_PressureVessel.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Axisymmetric Problems – Thick-Walled Pressure Vessel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/jpeg": "/9j/4AAQSkZJRgABAQAAAQABAAD/2wCEABALDBoYFhsaGRoeHRsfIiYiIiIfIiclJyAiMScyMC0nLSs1QFBCNThPOS0tRWFFS1NWW1xbNUFlbWRYbFBZW1cBERISGRUZLxsbL1c/NT1XV1dXV1dXV1dXV1dXV1dXV1dXV1dXV11XV1dXV1dXV1dXV1dXV1dXV1dXV1dXV1dXV//AABEIAWgB4AMBIgACEQEDEQH/xAAbAAEAAwEBAQEAAAAAAAAAAAAAAQMEAgYHBf/EAEIQAAEDAQUEBggGAQIHAQEBAAEAAhEDBBIhMVFBYXGRExYikqHSBTJTY4Gx4fAUFzNCUsHRI2IGFXJzgqLxQyQH/8QAFgEBAQEAAAAAAAAAAAAAAAAAAAEC/8QAHBEBAAMBAAMBAAAAAAAAAAAAAAECETESUWEh/9oADAMBAAIRAxEAPwD5+iIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAi9V1BtftKHef5U6g2v2lDvP8qDyqL1XUG1+0od5/lTqDa/aUO8/yoPKovVdQbX7Sh3n+VOoNr9pQ7z/ACoPKovVdQbX7Sh3n+VOoNr9pQ7z/Kg8qi9V1BtftKHef5U6g2v2lDvP8qDyqL1XUG1+0od5/lTqDa/aUO8/yoPKovVdQbX7Sh3n+VOoNr9pQ7z/ACoPKovVdQbX7Sh3n+VOoNr9pQ7z/Kg8qi9V1BtftKHef5U6g2v2lDvP8qDyqL1XUG1+0od5/lTqDa/aUO8/yoPKovVdQbX7Sh3n+VOoNr9pQ7z/ACoPKovVt/4AthMCpQ7z/Krfy4tvtbP3n+RB49F7D8uLb7Wz95/kT8uLb7Wz95/kQePRew/Li2+1s/ef5E/Li2+1s/ef5EHj0XsPy4tvtbP3n+RPy4tvtbP3n+RB49F7D8uLb7Wz95/kT8uLb7Wz95/kQePRew/Li2+1s/eqeRPy4tvtbP3n+RB49F6rqDa/aUO8/wAqdQbX7Sh3n+VB5VF6rqDa/aUO8/yp1BtftKHef5UHlUXquoNr9pQ7z/KnUG1+0od5/lQeVReq6g2v2lDvP8qdQbX7Sh3n+VB5VF6rqDa/aUO8/wAqdQbX7Sh3n+VB5VF6rqDa/aUO8/yp1BtftKHef5UHlUXquoNr9pQ7z/KnUG1+0od5/lQeVReq6g2v2lDvP8qdQbX7Sh3n+VB5VF6rqDa/aUO8/wAqdQbX7Sh3n+VB5VF6rqDa/aUO8/yp1BtftKHef5UHlUXquoNr9pQ7z/KnUG1+0od5/lQeVReq6g2v2lDvP8qdQbX7Sh3n+VB5VF6rqDa/aUO8/wAqdQbX7Sh3n+VB9DqAlpAzgxyWWzUiHtPRlsZkxpxM7OS/YdSYBJAUii3RZtXZifTUWyJh+WW1ZMERM8Nyjo6sEX2+Oa/UNNkxGa5qNY0SRh8T8IW9ZfnOZUklrhmc5OCXKn8h9/BfoUTTqA3cYN05iDoUJp37mF7TH7ncmo/Pc2oQ3tAHGfH6clz0VSB2xtnPHQr9V1No2KhtSmXXQMeBg6wcigxCnUEdoGN5x4+PPcpbSqA4vkYrcCzHLDPFReZeu9mYmJ2IMldryyGOuu1OPgu4Mt4H5LQSzD1ccsc1HSU9W80GYt7Xqn/qnD5/0uajakm65oG/VbGlhiLpnLFTLP8Abhv+9CpH4szrDcqfyH3ooiqZxaDOHBfodnQKJZE4RxVRhu1f5Nz5DTJSWVP5D5ff0W6G6DmoBYdgQYOjqfy8T/j4+Cksqzg5v3E/2t3Y0HNT2ZiBPFBgDKn8h9ngrSDdOsGFqhugXF+n/t5oKrE1wuh7rzsZIET8NistlGq57LhhozxIOa7ZUZsu89VZ041CxevlGNVnxnSqx5i6Yw2/D7+KqNKt/MffwVn4gRMiFPTYxhK0hRY8TecDllp/SqLa4yLTh46K01wMZEcVPS8EFdNlUHtPBGOzHb9Fw2lXAjpAcsT8M8Fd028YoK07RogqDK38m/fw+/BW0GvAN8g44QnTbwobXBmCMMCguVNnZUBffcHAns/7Rpv4/wCMZFbbIQVhEyI4oMdQGDBg7CRMfBVWRj202Co688ASd8ePFbnFs4gSfErm8z/brnsVRSiuvMicI1ldQ3QIM6K8lmOWGJx2ILhxgc0FCLQGt0Ch9xokgR8UFCK6mWOEtAI+KOfTDg03QTkJQUotbaTTsVFa00WVG0z67iMACYnInRNFaLTXNKm0vfDWjaUo9FUaHsIc05EHBPozIop26g55aA7B1ybpu3piJ4rVW6Km288ta3UmArMTHUi0W5LMi00ujeLzC1w1BkLhlag6o6k1zTUZF5oMls5SNimqpRWW2vSoNDqgwJgQCcYn+lNjq0q7L7BhMY4Yq5Obn4z5V8vHf1QWyZkjcF0rw+kXXcJmIxzUWupRosL6hDWjjidBvU1pitVn6QATsOuBMY79vNXhaqYpvaHtgtIkEHAhfgO/4os4rMaaNQUnuLW1sLpIMTEzG9TVx+nUZeETG8ZjgpaIAC2OpsAJIwAnCT4KRSborqOngEQclLQBkukUVWWNJBwkJVYHCCYyIIzBGRViIK6bboiSd5zKrNmbfvXj6wddwguAicp8VoXLskEOx2rM2yAOvSdwkQM/8nmrK1Qtuw28XGM4jAmfBV0LWHuu3SDE7olBy6wtLr22ZzGc5rgejKeOZnOXb5W1ZaltuueC31f9wk4A5bM0BthaDImc8xvn5p+BZ9lW0K18ExEGPD6q0zGGaDMLG2I/v71Kg2Fn2R97VH48SJa7HKB8/iFfQrB7bwkDLHBBV+Cb89o25o2xtGXKfj81ocYBOglZT6QaM2u27BsQT+AZ9kKTYmkyd+3XP5ruhab5waQInGMfgrkGUWIQQTM6YQNPAKfwTd+zbpkj7aGkyNsYYniRsGCtoVb7A6InZ8UFX4JmzDgcsZUmyg7T8CMOHIK6q+61zomATGsBZvx4yLXTEkDZu4oOvwbddkZjSFH4Fszpv3yumWwF9wNcDMGdmB/wtCDMLE37I+9qgWIYY5Y/HZ/a6rWkscG3SZjGdTA8VX/zFu1rvCeUoJbYWjbr4rttkaAY2iDliPsqKNsD3BsETlMaA/2tKDMbCzHfO3XNBYm/ZCmta7jiC1xA2jgT/SULWHm6AQc8YyQQLE3MTOGMjYpNjb8tughaEQZ/wbZBnEZYjDGVH4Fv2csIWlEGcWRoyMaYjDhyCj8Ez7O+VpRBl/Asx37wjrECSZzM+OPNakQZm2NoP13QgsTJn+98rSiCqnRDRAMDSQpfTDgQcuKsRBTSoBgIBwmcTMLl9ka54eTjhhIgxMSNua0IgluCyUfRlJld1cCajsydnAbFqXNRxDSQJgTEwg6qsa9pa7IpTY1jQ1uAGSpdaQMwZSnaQ4gQcTGKDGz0M1tQP6Q4PvxA/kXXZ0kr9C0UGVG3X4jHbGwj+yuK9dtNhe8w0DHaopVnOpXy0sME3TiRpK1a026zWsV4tpU2sENwEk/EmSuGWdrXlzTEmXAAdo6k5qt9quucC03WxLhvjZ8dVw23t2tcMY2aY7VlpprUg+MsJzEhU1KDm0ntYe0cohugw5ILWCGkAw5wGJg4xs25hRXtVx92B+2ZdBN4x2RGKsTieMPzfwFfO9Vnc9mB5r9J1kFWgynWJODb4B9YjYTotbV0tWvNupFYhmsdlZRpCkz1Rrmd5XnupzOkZ/qk0Wkm6QC44yG3tNi9Uixje4rqNDmluwgjYuhAXSIgihSgIoUoC5dkulBQcIsYp1haC7tlhIHrC6BGYC2jPn/SCEWWwUq7RV6d7Xy4ll1pEDmeSttgqGmei9eRGX9oLUUU70drOdmXwWSsyr014OdGwDEEaRsxkznvQbPvaiy+lqVV9nqtoG7VI7BmIdxT0cx7aZa8uMAGX+tJ9YTtxQak+9qo9IMquoPFBwbUg3TAOO6cJ4q2zteGMFRwc8ZkCAd8IOkWSmy0/inFxZ+Hui6Mb0y7ZriMdwwWp4MGM4McUEz94osFmo2roQKr2mrIlwwaewMv/JaGGtImIGeUnHdu8QgvSfvFY/TLK7rNUFmJFbC7BAOYnE7pVfoGnaG2YC1EmreMyQ7CcMQg/Q+9qfe1ZRSq9PJLw0EkQW3C27AaW5zOM+MYKfSDK7mN/DuDXXsS7K7Bw5wg0rkgEgmCRkdFzQFQUmCqQ6pdF4tEAujGAuX1KgJhs45YRGzFBbKSq7PUqEw9kCDj8coWhBXKclYodMGMDsnFByERl7C9E45ZBZ7WK0t6KYjGLox+IQaIU3SuKN662/623LP4K1BzdKXSu0QcXSl0rtEHF0pdK7RBxdKXSu0QcXSl0rtEHF0pdK7RBxB1SDqu0QVlmsKSyRB2rtEFZZOYBU3eC7RBxdKgs2wMF2iAApUIglFCIJRQpQZ+lduUdK7ciIHSu3KelduUIgdK7dyXVN7nAzGa4Kss+R4oOoO7kkHdyXaIKK94MJaJdGECcUpl5vXmgY4bxJz5Dmr0QfnmvUu3rrbvR3pw9fDDOVrxxwGemxWqHBB+f0toE9hh9c6TB7IzzIXVlr1nOIfSa1sYHU4LZdGvyS6NfkpjOfX55r1xLi1t1t8m8IJAJugRqIU9LaJ9RkSDF04C6SRM5zgt90a/JLo1+SrTFSr1ye1TaBjOuRyx1w3o6vWvnsC5LQOzjBPaOa23Rr8kujX5IPz/AMVXmOiGTjl/1XYx3DnsQWutfg0gGktDSczIlxidi/QujX5JdGvyQY7TWqtcOjYHC6dn7pEYzpPJU0rXaHf/AItwcAf72+K/SujX5JdGvyUxmY+sdGpWLm32sALJMD1XYYTKrfWtAeAKbHNLwJ/i28QScdIK/QujX5JdGvyVWH59K0V2i9UY2AyXXWGb2MgCTOQ5rTUqOBbDZBmYGWGCvujX5JdGvyRX5jbTarpcaLJu4Mj916MTP8YMR8VroVKhvX2NbBw3j7/ztWi6Nfkl0a/JBxJ0HJJOg5Lu6Nfkl0a/JBxJ3ckk7uS7ujX5JdGvyQcXju5Kb7ty6ujX5JdGvyQc33bkvu3Lq6Nfkl0a/JBzfduS+7curo1+SXRr8kHN925L7ty6ujX5JdGvyQc33bkvu3Lq6Nfkl0a/JBzfduS+7curo1+SXRr8kHN925L7ty6ujX5JdGvyQc33bkvu3Lq6Nfkl0a/JBzfduS+7curo1+SXRr8kHN925L7ty6ujX5JdGvyQch7tyl7nAE4LoNGqVfVPBBQKrt3JT0rty5bkpQT0rtydK7coRBPSu3Lk2hwc1sHGcQ0kCNTsUroIIREQEREEFWWfI8VWVZZsjxQXIiIKbVVLKbnCMBtySlWvXuyRdMY7cSJHJWpCDF+LfGQ/S6Sdk6blqLpncY+Snom3bt0XdIEclLhIQfmt9IyTLIEvgYS660HbxI+CD0m05CILJvf7ruG4w7boV+jdKgM4Y+KDDT9Ih0gMeCC8CRndEk+I5rmp6RLHQ6mS2QAW4zlJ8St/R4zAmInbGikMjKEGAekHYk0iIDpgzBF7A8Q3xWllQkkXciBOw5TG1X3SoDIwEIMj7ZdzYSRMxlInD4x4hcst979t3Fox0d/a23Sl0qNbHpkbbC5he2mYEiDnOzAbFaK3YY6MDE7pH+YVrsASSABmVN0qpsMA9Iy0uDDduyCYzxMeCvFo7UBs9oDCMBdBk81oDOCXOChselAr/wCnfLSYkQ0STBjAKXVPUgesYx2YE88FaRdGwAfABDTnMA7fiqKqVWWhzoAgzxBiRuVD7aWvcHU3ReutiMcBJz35aA6LY2nAAEADIDYpulEfnVfSlwG9RqYCf2mdBgd6tp2wuMXLvaIxOJAbeHxK2XSlzgg/P/5mLoIpvJLQ6BGRMf0eStq2wN6M7HTIjEQwu/oD4rUKeJMCTmdp4qbpQfmf827Jf0L4uhwGF4y4tOGUCM52hXutwF43SA2/MxiGtns8f6K23SuXU5zAPHFBgPpIQeyWxBl2UXgNnHA7YOin/mU9JdpPJZfxMAG7GPxn5r9C6UulB+efSOjCYAMCJg3pP/qOaN9JTfIpVCGFwOAE3Yyk7Z8Cv0LpUXSg/Pb6UBH6NThAkYTru+KsFvGPZkjpMojsbOJWy6VN0oPz3ekIeWBhcWkBxbEepekeGeasp2u8MLuJqN4FswTuw8QtTacZADhghpzmAcIx00QZvxR/jspmNsPMRxCupvvA7iRyK76PcNfjqpDI0QQim6dyXTuQQim6dyXTuQQim6dyXTuQQim6dyXTuQQMwuqvqnggbilX1TwQZ25KVDcl0ghERAXQXK6CCEREBERBBVlmyPFVlWWbI8UFyIiAiIgKFKhBD3BokkAalSsLPR4v1L2LSZbBIImCQTxGHFUss9o6VzyR6xETg4BpuncAT9xiH6iLJSoVAHA1CTdaLxiS/GTAEbQMtisp06gYQ6pLpwdAMDhAQVUbYatT/TANIGHOnMlsiBpi3mti/K6VtNpu1mgXohrGiTN2Y28dy6fZ65rOcXC4C0AD9zb0kAbMI5b8A/TXFUw1xxy/aJPwGqz0KNUOBc+TddewEXiRAERgO1zXdGnUDCHPF+TDg3IToTogy1BUc2kCWkbQ9xY55GMEXdkfFK9rrE02sYJOOEkECQRMCP25jarKvo9zyC6oCRH7BsM5Sumh14sFdt5oEtDWyBswQcUq1UhpIESwZYuJALzuAk/ELes4o1fa/wDoFaWuvYOF3S7jzn+kGS2h5ewdi5ibpeWlxjcNiqtNrr3qbadMSRf2kOGIIJIERLT8VdXsBe++XjdLZgAyBnqodfaQ0VhMgQGtwndKCaVaqS2R2bzQMMXC7LnHTH5b1tWdtGrImrI/6BipLKnSg3h0cEXYMzhtnigptT3dKwAgNBmLwBfpyKoZbaprOFwimHXDh6pEmd96Wjd8VoqWEuqdJfE7OyDGi5Dnl4aKwOeIY0gERgd+KCyyVajj/qNDTca4gYhrjPZnbsWpVUqbwe0+8NLoCtQEREBERAREQEREBERAUqFKAiIgIiICIiAiIgLir6p4LtcVfVPBBnbkuly3JdICIiApChSEEIiICIiCCrLNkeKrKss2R4oLkREBFxVqBjS50wM4BPgEbUaZggxgY2HQoO0VH4tkTJ9W96pwbqrC7Q7YQdos4rtJAvtkkgCRJIzC7vYxOKC1V1qQexzDk4FpjQiFX+IbBPSNgZm8MPuFD7S0BpLsHEBpHaknggmz2NlMktGJEEnExJMc3E/FXrP+JZ/MbZMiBETJ2ZhPxLcIeDJA7OOeUxkDqg0IszrXTAk1G7f3DZmn4tkgXxiQ0YjEkSBPBBpXLKYaXEDFxk7zAHyAVJtTPas2/ubszQWlntG7B6wzOQQaEWf8UyJviIJzzAzIG1DamBwaagvExEiZ04oNCqZZ2tLnAYudeJOsAfIBcsrh0Q6ZEjeJjBci0tLi0OxBg7jj/goNKKuTqpk6oO1XQoNptDWiI5naSUk6pJ1QWKVVJ1UydUHaKuTqknVBYirk6pJ1QWIq5OqSdUFiKuTqknVBYirk6pJ1QWKVVJ1STqgtRVydUk6oLEVcnVJOqCxFXJ1STqgsRVydUk6oLFxV9U8FAJlTV9U8EGduS6XLcl0gIiIIXQXK6CCEREBERBBVlmyPFVlWWbI8UFyIiDirTD2lpmDnBI8QoZRa2YETnzJ/srtEGb8E27GN65cvbuGSvIjmF0ocg/LDLMwhzTduvc7C/if3YbVe600ms6cYh20Ay6NkfAq38IyZuDbpEmJMa4BDZWEQabSMMIGyY+Z5oMTqdAnCXfqdntYu/dicsiBiBmuy+ztawE3QHNLRLsHESPBx+wtX4Nns27dg2580/CsgDo2wIjAYRlyQYTToNll9zb96BJJuw0GMz+wR4K6jZ6LvVDh6jh2njDMHPUYjbGK0mysObAc84OBzHBTSoBghrYAAHwGQRFNOwUmtLWtgHPtOnnMroWKmNh2ZucQIECMdFfB0UwdEVkZ6OotaWhmBicSZgznOpRtgpDEN8XYcMcFqg6KYOiDMLEwAgXsQRN5xzJO06ko2xMBLiCXEzJJ1nAbMVpg6KIOiDP8AhGgQ1zhg4A3iSLxkkEyZ+S6FlYHXgIJdeO8iY8TOCvg6JB0QQimDokHRBCKYOiQdEEIpg6JB0QQimDokHRBCKYOiQdEEIkHRTB0QQiQdFMHRBCKYOiQdEEIpg6JB0QQimDokHRBCKYOiQdEEIpg6JB0QQimDokHRAGYU1fVPBQAZU1fVPBBnbkuly3JdICIiAs9S0lr7sDYtCwWj9Xkg3oiICIiCCrLNkeKrKss+R4oLkREFVoa4thueG0iRpIyVDadaW3nDA4wYvCDnhwWtFMXWDoql3bHRxEm8H/8A3atrZxn4LpEiCZ1+W211jVd2CKYddIjERJEa3pA3QtNGpWg3mi9cBhuQeZ7IJz2LWiqMVWpWNKbpa6+0GIJDJF4gScYlW0qh7TQCSACLxjOYGWGS0IgwuqVjSpktuuLhfu4kMjEgRnMDgpFZ5d0bADB7TnGbrSCYIEY5BbUQZHipTYxtNrXQQDm0ATpioe6p0r5vBgDSyIgnEuvHkOC2JKDFTv1iHOaGsBa5oIlx7IM7sSR8FTUtdamHEsloc7E/xvOiP/XDbK/TRBhFeqbt1olzWHI3QTN864QOa66Wr0YLWBxvkYmOxfgHI/txWyUlBjqOf0rr14MDWltwTJk3pPLD7Fb7TVN19NhLXUwYIwBMn/A+Mr9BEH59K2PeCYECpdMTg3HbtxAHxV1lqVCR0gANxpIGTXEnsztwWpEBESUBElJQEREBESUBElJQESUQEREBESUBElJQESUlARJSUBElJQFzV9U8F0uavqnggztyUqG5KUBERAWG0fq8luWG0fq8kG9ERAREQQVZZsjxVZVlmyPFBciIgptNW4xzsMBtMDmlKvevYEXTGMY4kSOStSEGQ2l929db+nfOJw3ZceS0F0zuMfJOhbdu3RdyjYpfkgyWW1NqThdN4gA4FwGMidxUMtQNWpTIDQy7DpzJjCNmJw1Wm6JBgSMjtCFoOYB2466oKalYCo1nZAIJkmMtgG3euLXaeja1wAcHHWJF0nDUmMBvWkgEEHEHMHapCCtzwHBsHHbGA4lUvtjA9zbpkFowAzcDHyhalDmgwSAYxG4oMlK3sdTc+65t2TdddDiBGIxjaFy/0iwR2Hu9XIDCRIOJ0W5Q1oEwAJMnedUH549KNLJuQ64XYlt2QCcwcsI3bcVe21tL7rRjeaJwiC0ukRuBWpQROaCh1oDXhrgBgSXbBnAngCUtNe4WC7N5wbnlv3q8gEQRI0KlF2FNWqGsLhBGAzwziSdAui8XSRBIGIB2xMLsNAEAQNyBoGQzz3oMzbQTTNS62ACc8DGhjJaYQgHPEKUSUQl0KUQRdCQFKIIhIUogiEhSiCISFKIIhIUogiEhSiCIS6FKIIuhLoUogi6EuhSiCLoS6FKIIuhLoUogNzC6q+qeC5GYXVX1TwQZ25LpctyUoJRQiCVgtH6vJblhtH6vJBvREQEREEFWWbI8VWVZZsjxQXIiIIRFKAuX5cvmuly/Lkg/NpvtJ7RDQO1hGO4+tnPPctINTo2YC/AvB2zDHJZm2euHEl4Ml/7spaA2Bd1E7p2qBRtB9YtMFkdqAYulxENn+Q3zs2hfaBVL2hkBkEOMiQSMCBtgxzXBfXabsBwljQ6MTI7RMHYcdmi5p07RJDngg3yCIwH7W65kn4BQ+zV5vU6l2SOy4yGjDAZ6Jpqyk6sab+kaL/au9HA2CMznM8lTWtNobd/0xBc0ereIEdo4O/x8V0LLWExVMw4CTMzegnDCJHLatVOmdpMYXW/xAAww/wAoKgasGTAh04S6ZMXd2UYGVaS+62A29hek4DWFQ+jVJuh0MAznHM7s8uW9dNpVAe0QRLTgcTAg4QIxgprWIrsrGr2HXadyJhph+OOv8V2KdUnF4ADhENxc0ZyZ28App03gtvOLiAbxwxOGEc1Yxrg5xLpafVERHx2pqM5pVvV6Q+q7tw2AZwwzUOqVw66GggkBpjMXZJMHDEbswuqVN0gOn1XhxnMl2EY4bV3aKLnlsOuhrpwPrC6cDhrHJNJVTXNF4c0dLBulkATszP3CrrWiu17R0YuueBIaSQ3bMHSMcNui66C0Fp/1odAAi6RN4ySbukbP8qKdntF8F9UEBzTDcJbdIdPZ1gx4hEdg1ujJOBuZR2jUxnbETGzLatLy+W3QIntSchGz4rGWV2tF0iQwCCQZqXsXcDjt+CuoU6rXOvvD2k4ThdGOg4eKBVZUc5wDi0FoukXSAdsg4z4KQyqSJeGgOOAbm3ECTPA7Fl/DWi6B0jWgMAhrjN8On1i3KMMtuSup0q4fJqNLb5MR+zYMs/vbhdFlFtQOAc4uAZi6GgOdOgxC6d0l2pgA7G5B3YTO2Vk/D1zAv3GhrdoJLhEnLcdu3JXU6dS92oPbLpDp7FyIIgbdmPFQXWYPuC/N7HOJiTExhMRkrVkbRfd2zcAz/dOB+G3Va0BERAREQEREBERAREQEREBERAREQEREAZhdVfVPBcjMLqr6p4IM7clKhuS6QQilEELDaP1eS3rBaP1eSDeiIgIiIIKss2R4qsqyz5HiguRQiAiIglQQiIKOnp3yye0IwnWY+RXTKjHTDgYxkHCOPwKzN9FMDi+TevFzT/CRBA4q9ljYAQBDbgYAMIaJwEcUFgLSJnDiqmWmk4NIdg9pe04gFuGPiFLbIwMLBN0mT2iZ3Y7FwPR7AGAT2Ju45SQf6CC01GY9oYRkZOOA5o17HGAQTuOmf9c1moejgwOE4S0tAHq3XSBqf/uqup2NjbsDBrS0A6OILp1JgIOqtVjKbqhPZaCSRJyzVItlNzQWEEl12HEtODrrsOIPFdfgKfRGkBdaTJiMcZXdmsjKU3BnmSZP/wAzw3oJbVYZF6CCRiYxGaGrT/kIiZnCJiZVFf0ax18iQ5+04icYw3XjzXQsLYbeN6AJwAvQ68JjQoLw9kTOExntmI5ripWY113Eui9DQSQ3VcvsTHNLXC8C6/jrMrv8M2/fEtN0NMZFoJIHwk80HLazXVLrS0iCSQ6SDMRC7c9ga509lsyZyjP5KaNFrGBjRDQIVJ9H0+iqUgLral4Ou4Eh0z80HdnqsqA3ZwMEGQQYB+RCuuBcUaIZejNxvOJzJiJ5AKxBFwJcClEEXAlwKUQRcCXQpUoOboS4FKIIuD7KXB9lSiCLgS4FKIIuBLgUogi4EuBSiCLgS4FKlBzcCXAukQc3AlwLpEHNwJcC6RByGhRV9U8F0uavqnggztyXS5bkukBERAWC0fq8lvWC0fq8kG9ERAREQQVZZ8jxVZVlnyPFBRYnl7q14k3ahAxyEZLVcGp5lY/R3r1/+6fkFuUrxmvHJYNTzKXBqeZXFpo9IxzcMRGIvDklKiW3u0TeM4zhicvDkq06hsTOGt4o5sa8yqTRPQdHDb1y7nhMRoriIB3n/CDn7zKRx5lYKtjqQbtXtHpbsh37hgAb2EHcfguGejKl0f8A9D2uhswXOGEZSdyD9KOPMpHHmVnpUGhzT0jnEThfOJmCYn4RsVzWQ4uvEzsOQ4IOvvMpHHmV+fVsdSCG1u0W1AJDp7RGRnCPj8FDPRtS63/XeCLpIBc4YbASUH6MceZT7zKzVKZdJe7ozdwDXmANrtmvwXD7K4uhlSAHtccXF10Mi7nqAcZzOCDX8fEqY48yvz22UMaCKgLBScIdJBBDcZnAYabV+g0QAJmBE6oEceZSOPMqUQRHHmUjjzKlEERx5lI48ypRBEceZSOPMqUQRHHmUjjzKlEERx5lI48ypRBEceZSOPMqUQRHHmUjjzKlEERx5lI48ypRBEceZSOPMqUQRHHmUjjzKlEERx5lI48ypRBEceZSOPMqUQRHHmUjjzKlEERx5lI48ypRBAHHmVLj/p/D+kGaH9P/AMf6QUtyXS5bkukBERAWC0fq8lvWC0fq8kG9ERAREQQVZZ8jxVZVlnyPFBm9HevX/wC6fkFvWD0d69f/ALp+QW9ZrxmvEIiLTQoflyXS5flyQfk1LNSIc1tWCemx7JAkAP2ZDDaqKVgshaP9Vsw2SxzWziIw5L9CrYKbmuHaE3snOzcIJiYPxXI9F0boBbeiMSTJjgg4oU7OKrbjm35qEAXdpBcMBrG9a2MaHuI9YxOPLBBRaC0gRdkCCcjmN/xXQYAS4AScztKo/NqWaiWlrasS2rj2SACRe2bDG1VUrFZbrYqtkXJLHNaCf24DwW+r6OpuaW9oSHDBzv3Z4EwfioHo2lDQW3i2MSTMjaVBVa2UsTaHUwbjgYAEsOeBk7NiipQpOLZqAND2EYNEu6OAAeF04LY2hgb5LyRBJgdnTCNSuHWNhIJvYEOi8cSG3cdpwjPRBlo06IaDSf8A/m+AAD2ezJiJmQMNZX6DAAABlAjgqjZGZiWm6WiHHAEAZZbByVwAAgZBBKIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiIAzQ/p/wDj/SDND+n/AOP9IKW5LpctyXSAiIgLBaP1eS3rBaP1eSDeiIgIiIIKss+R4qsqyz5HigzejvXr/wDdPyC3rDYmlrq14ETUJGBxEZrXfG/kVmvGa8dIub438inSDfyK007XL8uSi+N/IqHPG/kUGBvpNpJEQJeBLmybrQcp3nhGMKB6VYcoGLJvEAdqN+fa25wdFun7gqBAnDPE4HEoMtP0kx0w12BeMs7ok/1zXFT0lcdD2GJABbjeOGXNbYEzGMRN0zGikYCBlwKIwt9JZk0yIDpxaYLb2BjW6tTahJIu5EAnYcpjbtVk/cFQMBAGA3H/AAi6zvtobgR2occxEtXB9JMkgNcYiSIjETOa2T9wVy1jQSQIJzIaZPgo1E1VMtYdSdUaDDZwO2FP4g3ouk4huH/TJOOxduptJktBMRJactMlLWgZCMth2ZbETYU/ijAIbEsc8TujOOKtfVi7g43iBgJjjopcwHMThGLTlpkup+4KpsKq1YtvQJusLsZ2bPmu7+sCT2dThMY7c0exrsxOzFpy0XU/cFE2GKn6RkNvU3B7gHBoLZgmBt355KP+atv3Cx4deuCQMTjv3LdP3BSePIoMTfSI6PpC2AGBxEtJBLiImYwgrp3pBoeGXHkl10kAQ04Z46EH4rVtmMdYP+EaAMhGJOAOZzOSDJaPSLKd4khwDA8XSC4gujLTLFQPSbbwaWOBLnNGLdgmTjhIgwVtn7gpPHkUGQ+kG3C+6boDTBgO7TiMjllx2Lip6VY112C43nthpaYugEzjnjktu2duUwVy6m0kOLQS3IlpkcMEGSn6Va43W03l0Ax2Qc4jOFe20k3YaQHOc0Xo2AmcNnZKvnjyKggHMbsjkgobapAMZtY/PUwR/jVWWZ5dSY45ua0niQuro00/admWxSwBoAAgAQBBwCI6RRP3BSfuCipRRP3BSfuCglFE/cFJ+4KCUUT9wUn7goJRRP3BSfuCglFE/cFJ+4KCUUT9wUn7goJGaH9P/wAf6UA8eRUu/T+H9IKW5LpctyXSAiIgLBaP1eS3rBaP1eSDeiIgIiIIXVOoGgyDmoUILOnGh8E6caHwVaILOnGh8E6caHwVaILOnGh8E6caHwVaILOnGh8E6caHwVaILOnGh8E6caHwVaILOnGh8E6caHwVaILOnGh8E6caHwVaILOnGh8E6caHwVaILOnGh8E6caHwVaILOnGh8E6caHwVaILOnGh8E6caHwVaILOnGh8E6caHwVaILOnGh8E6caHwVaILOnGh8E6caHwVaILOnGh8E6caHwVaILOnGh8E6caHwVaILOnGh8E6caHwVaILOnGh8E6caHwVaILOnGh8E6caHwVaILOnGh8E6caHwVaILOnGh8E6caHwVaILOnGh8E6caHwVaILOnGh8E6caHwVaILOnGh8FDqwIIgrhEAKURAREQFS+zBzr06K5SEEIvmHX23e67n1Tr7bvddz6oPpylfMm/wD+gW0D1aB3ljv6cp/MG2fwod13mQfTEXzP8wbZ/Ch3XeZPzBtn8KHdd5kH0tF80/MG2fwod13mT8wbZ/Ch3XeZB9MRfM/zBtn8KHdd5k/MG2fwod13mQfTEXzP8wbZ/Ch3XeZPzBtn8KHdd5kH0xF8z/MG2fwod13mT8wbZ/Ch3XeZB9MRfM/zBtn8KHdd5k/MG2fwod13mQfTEXzP8wbZ/Ch3XeZPzBtn8KHdd5kH0tSvmLv+P7achRHBh/srnr7bvddz6oPqCL5f19t3uu59U6+273Xc+qD6gi+X9fbd7rufVOvtu913Pqg+oKF8w6+273Xc+qkf8fW73Pc+qD6ei+YdfrbpR7h/ynX626Ue4f8AKD6ci+Ydfbd7rufVOvtu913Pqg+novmHX23e67n1Tr7bvddz6oPp6L5h19t3uu59U6+273Xc+qD6ei+Ydfbd7rufVOvtu913Pqg+novmHX23e67n1Tr7bvddz6oPp6L5h19t3uu59U6+273Xc+qD6ei+Ydfbd7rufVOvtu913Pqg+novmHX23e67n1Tr7bvddz6oPp6L5h19t3uu59U6+273Xc+qD6ei+Ydfbd7rufVOvtu913Pqg+novmHX23e67n1Tr7bvddz6oPp6L5h19t3uu59U6+273Xc+qD6ei+Ydfbd7rufVOvtu913Pqg+novmHX23e67n1Tr7bvddz6oPp6L5h19t3uu59U6+273Xc+qD6eoNVoc1pOLpgQcYzXzHr7bvddz6qR/x9bvddz6oPLoiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIg//Z\n",
      "text/html": [
       "\n",
       "        <iframe\n",
       "            width=\"400\"\n",
       "            height=\"300\"\n",
       "            src=\"https://www.youtube.com/embed/Z2eaaBGDOy4\"\n",
       "            frameborder=\"0\"\n",
       "            allowfullscreen\n",
       "        ></iframe>\n",
       "        "
      ],
      "text/plain": [
       "<IPython.lib.display.YouTubeVideo at 0x298aa7288b0>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.display import YouTubeVideo\n",
    "\n",
    "YouTubeVideo('Z2eaaBGDOy4')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will follow the development of the ideas presented in Timoshenko and Goodier.\n",
    "\n",
    "Earlier, we had considered various definitions leading to the biharmonic equation and stresses in polar coordinates. \n",
    "\n",
    "In this file, we are going to consider the set-up required for axisymmetric problems with particular focus on thick-walled pressure vessels. \n",
    "\n",
    "Before proceeding, however, note that in this file, we do not go into writing all those definitions again. Instead, we utilize those definitions packaged in a separate python file `polarUtilites.py`; this file is stored in the same folder as the Jupyter Notebook files. So, before proceeding, in order to make those definitions available for the current problem solution, we first need to import that file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from polarUtilities import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we will consider a $\\theta$-independent Airy stress function, i.e. $\\phi(r)$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\phi{\\left(r \\right)}$"
      ],
      "text/plain": [
       "phi(r)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "phi = sym.Function('phi')(r)\n",
    "display(phi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's quickly check what the polar biharmonic of $\\phi(r)$ looks like now that we dropped the $\\theta$ dependence:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{d^{4}}{d r^{4}} \\phi{\\left(r \\right)} + \\frac{2 \\frac{d^{3}}{d r^{3}} \\phi{\\left(r \\right)}}{r} - \\frac{\\frac{d^{2}}{d r^{2}} \\phi{\\left(r \\right)}}{r^{2}} + \\frac{\\frac{d}{d r} \\phi{\\left(r \\right)}}{r^{3}}$"
      ],
      "text/plain": [
       "Derivative(phi(r), (r, 4)) + 2*Derivative(phi(r), (r, 3))/r - Derivative(phi(r), (r, 2))/r**2 + Derivative(phi(r), r)/r**3"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "polarbiharmonic(phi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's display this in a better fashion:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\nabla^4 \\phi(r) = \\frac{d^{4}}{d r^{4}} \\phi{\\left(r \\right)} + \\frac{2 \\frac{d^{3}}{d r^{3}} \\phi{\\left(r \\right)}}{r} - \\frac{\\frac{d^{2}}{d r^{2}} \\phi{\\left(r \\right)}}{r^{2}} + \\frac{\\frac{d}{d r} \\phi{\\left(r \\right)}}{r^{3}}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from IPython.display import Math, Latex\n",
    "display(Math(r'\\nabla^4 \\phi(r) = {}'.format(sym.latex(polarbiharmonic(phi)))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we would like to find the solution of the biharmonic equation, $\\nabla^4 \\phi = 0$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\phi{\\left(r \\right)} = C_{1} + C_{2} r^{2} + C_{3} r^{2} \\log{\\left(r \\right)} + C_{4} \\log{\\left(r \\right)}$"
      ],
      "text/plain": [
       "Eq(phi(r), C1 + C2*r**2 + C3*r**2*log(r) + C4*log(r))"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eq = sym.Eq(polarbiharmonic(phi),0)\n",
    "sym.dsolve(eq)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That's it! That simple one line gives us the solution to the fourth order equation. Anyone who has tried to solve the fourth order equation by hand can attest to the fact that how much labour this step is saving us. I had emphasized this point heavily during my lecture so that the students appreciate the power of SymPy.  \n",
    "\n",
    "Anyway, now that we have the solution, we will rewrite it in the notation used by Timoshenko and Goodier. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\phi(r) = A \\log{\\left(r \\right)} + B r^{2} \\log{\\left(r \\right)} + C r^{2} + D$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "A, B, C, D = sym.symbols('A, B, C, D')\n",
    "\n",
    "phi = A*sym.log(r) + B*r**2*sym.log(r) + C*r**2 + D\n",
    "display(Math(r'\\phi(r) = {}'.format(sym.latex(phi))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To determine the corresponding stress fields, we use our earlier definitions. We just have to call them:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\sigma_{rr} = \\frac{A}{r^{2}} + 2 B \\log{\\left(r \\right)} + B + 2 C$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\sigma_{\\theta\\theta} = - \\frac{A}{r^{2}} + 2 B \\log{\\left(r \\right)} + 3 B + 2 C$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\sigma_{r\\theta} = 0$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sigmarr = sigma_rr(phi)\n",
    "sigmatt = sigma_tt(phi)\n",
    "sigmart = sigma_rt(phi)\n",
    "\n",
    "display(Math(r'\\sigma_{{rr}} = {}'.format(sym.latex(sigmarr))))\n",
    "display(Math(r'\\sigma_{{\\theta\\theta}} = {}'.format(sym.latex(sigmatt))))\n",
    "display(Math(r'\\sigma_{{r\\theta}} = {}'.format(sym.latex(sigmart))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we are investigating the stress field for a geometry where there is no hole around the origin then we must set $A=0$ and $B=0$. Otherwise $\\sigma_{rr}$ and $\\sigma_{\\theta\\theta}$ will become infinite. In such a case, we have the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\sigma_{rr}^{\\rm no-hole} = 2 C$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\sigma_{\\theta\\theta}^{\\rm no-hole} = 2 C$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sigmarr_nohole = sigmarr.subs([(A,0),(B,0)])\n",
    "sigmatt_nohole = sigmatt.subs([(A,0),(B,0)])\n",
    "\n",
    "display(Math(r'\\sigma_{{rr}}^{{\\rm no-hole}} = {}'.format(sym.latex(sigmarr_nohole))))\n",
    "display(Math(r'\\sigma_{{\\theta\\theta}}^{{\\rm no-hole}} = {}'.format(sym.latex(sigmatt_nohole))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above formulae show that the radial and hoop stresses correspond to a state of uniform tension or uniform compression. \n",
    "\n",
    "Next, if there is a hole then we can have various other solutions. Suppose we take $A \\neq 0$ and $B=0$. Then, we have the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\sigma_{rr}^{\\rm hole} = \\frac{A}{r^{2}} + 2 C$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\sigma_{\\theta\\theta}^{\\rm hole} = - \\frac{A}{r^{2}} + 2 C$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sigmarr_hole = sigmarr.subs([(B,0)])\n",
    "sigmatt_hole = sigmatt.subs([(B,0)])\n",
    "\n",
    "display(Math(r'\\sigma_{{rr}}^{{\\rm hole}} = {}'.format(sym.latex(sigmarr_hole))))\n",
    "display(Math(r'\\sigma_{{\\theta\\theta}}^{{\\rm hole}} = {}'.format(sym.latex(sigmatt_hole))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Students who have done an undergraduate course on Strength of Materials/Mechanics of Materials/Mechanics of Solids can identify that this form of the stress fields emerges in the problem of a thick-walled pressure vessel. \n",
    "\n",
    "So we will now focus on the thick-walled pressure vessel problem. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solution to the thick-walled pressure vessel problem\n",
    "\n",
    "Let the inner radius be $a$ and the outer radius be $b$. Then the boundary conditions are:\n",
    "\n",
    "\\begin{align}\n",
    "\\sigma_{rr} &= -p_{\\rm in} \\quad \\text{at $r=a$} \\\\\n",
    "\\sigma_{rr} &= -p_{\\rm out} \\quad \\text{at $r=b$}\n",
    "\\end{align}\n",
    "where $p_{\\rm in}$ and $p_{out}$ are the pressures inside and outside the pressure vessel, respectively. We will need to solve for the two unknown constants $A$ and $C$ using these two boundary conditions. We can again implement this solution in SymPy. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle p_\\mathrm{in}$"
      ],
      "text/plain": [
       "p_\\mathrm{in}"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle p_\\mathrm{out}$"
      ],
      "text/plain": [
       "p_\\mathrm{out}"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pin, pout = sym.symbols(r'p_\\mathrm{in}, p_\\mathrm{out}')\n",
    "display(pin, pout)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "a, b = sym.symbols('a, b', positive=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We set up the two equations from the boundary conditions as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{A}{a^{2}} + 2 C = - p_\\mathrm{in}$"
      ],
      "text/plain": [
       "Eq(A/a**2 + 2*C, -p_\\mathrm{in})"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{A}{b^{2}} + 2 C = - p_\\mathrm{out}$"
      ],
      "text/plain": [
       "Eq(A/b**2 + 2*C, -p_\\mathrm{out})"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "lhs1 = sigmarr_hole.subs([(r,a)])\n",
    "rhs1 = -pin\n",
    "eq1 = sym.Eq(lhs1,rhs1)\n",
    "display(eq1)\n",
    "\n",
    "lhs2 = sigmarr_hole.subs([(r,b)])\n",
    "rhs2 = -pout\n",
    "eq2 = sym.Eq(lhs2,rhs2)\n",
    "display(eq2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we solve the above two equations to obtain the values of $A$ and $C$ as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle A = \\frac{a^{2} b^{2} \\left(p_\\mathrm{in} - p_\\mathrm{out}\\right)}{a^{2} - b^{2}}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle C = \\frac{- a^{2} p_\\mathrm{in} + b^{2} p_\\mathrm{out}}{2 \\left(a^{2} - b^{2}\\right)}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "soln, = sym.linsolve([eq1,eq2],[A,C])\n",
    "\n",
    "Asoln = soln[0]\n",
    "Csoln = soln[1]\n",
    "\n",
    "display(Math(r'A = {}'.format(sym.latex(Asoln))))\n",
    "display(Math(r'C = {}'.format(sym.latex(Csoln))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To obtain the solutions to the stress field, we just have to make the appropriate substitutions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\sigma_{rr}^{\\rm hole-soln} = \\frac{a^{2} b^{2} \\left(p_\\mathrm{in} - p_\\mathrm{out}\\right)}{r^{2} \\left(a^{2} - b^{2}\\right)} + \\frac{- a^{2} p_\\mathrm{in} + b^{2} p_\\mathrm{out}}{a^{2} - b^{2}}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\sigma_{\\theta\\theta}^{\\rm hole-soln} = \\frac{a^{2} b^{2} \\left(p_\\mathrm{in} - p_\\mathrm{out}\\right)}{r^{2} \\left(a^{2} - b^{2}\\right)} + \\frac{- a^{2} p_\\mathrm{in} + b^{2} p_\\mathrm{out}}{a^{2} - b^{2}}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sigmarr_hole_soln = sigmarr_hole.subs([(A,Asoln),(C,Csoln)])\n",
    "sigmatt_hole_soln = sigmarr_hole.subs([(A,Asoln),(C,Csoln)])\n",
    "\n",
    "display(Math(r'\\sigma_{{rr}}^{{\\rm hole-soln}} = {}'.format(sym.latex(sigmarr_hole_soln))))\n",
    "display(Math(r'\\sigma_{{\\theta\\theta}}^{{\\rm hole-soln}} = {}'.format(sym.latex(sigmatt_hole_soln))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These solutions match with those given in Timoshenko and Goodier!"
   ]
  }
 ],
 "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}