CurvedBar_PureBending.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Pure Bending of a Curved Bar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/jpeg": "/9j/4AAQSkZJRgABAQAAAQABAAD/2wCEABALDBgYFhwaGRoeHRwfIiclHyIhIyclJyctLyoxMC0oLSs2PVBCNThLOS0tRGFFS1NWW11bNUFlbWRYbFBZW1cBERISGRYZMBsbMGRANT1XV1dXV1dXV1dXV1dXYF1XV15XV1dgV1dXV1dXV1dXV1dXV1dXXVdXV1dXV1dXV1dXV//AABEIAWgB4AMBIgACEQEDEQH/xAAbAAEAAgMBAQAAAAAAAAAAAAAAAwYBAgQHBf/EAEgQAAIBAgIFBQwHBwMEAwAAAAABAgMRBCESMUFRkhdTYXHSBRMVFiIyNFSBkZOxBhRScqGywQcjJEJi0fAzguFDRIPCc6Kz/8QAGAEBAQEBAQAAAAAAAAAAAAAAAAECAwT/xAAlEQEAAQIGAwEBAAMAAAAAAAAAAREhAhMxQVGhAxLhYZEiMnH/2gAMAwEAAhEDEQA/APPwAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAWyX7PcYtdbDLrnPsm3J3jedw/FPsAVEFu5O8bzuH4p9gcneN53D8U+wBUQW7k7xvO4fin2Byd43ncPxT7AFRBbX+zzG87h+KfYM8neN53D8U+wBUQW7k7xvO4fin2Byd43ncPxT7AFRBbuTvG87h+KfYHJ3jedw/FPsAVEFu5O8bzuH4p9gxyeY3ncPxT7AFSBbuTvG87h+KfYHJ3jedw/FPsAVEFt5O8bzuH4p9gzyd43ncPxT7AFRBbuTvG87h+KfYMcnmN53D8U+wBUgWmp9A8VDzq2GXXKfYNqf0Axc1eNbDNdEp9gCqAtvJ5jdXfcPf70+wH+z3GLXVw3HPsAVIFtf7PcYtdXDLrnPsDk8xvO4fin2AKkC28nmN5zD8U+wOTzG85h+KfYAqQLbyeY3nMPxT7A5PMbzmH4p9gCpAtvJ5jecw/FPsDk8xvOYfin2AKkC28nmN5zD8U+wOTzG85h+KfYAqQLav2d43ncPxT7Bnk7xvO4fin2AKiC3cneN53D8U+wOTvG87h+KfYAqILdyd43ncPxT7A5O8bzuH4p9gCogt3J3jedw/FPsDk7xvO4fin2AKiC3cneN53D8U+wOTvG87h+KfYAqILdyd43ncPxT7A5O8bzuH4p9gCogt3J3jedw/FPsDk7xvO4fin2AKiC3cneN53D8U+wOTvG87h+KfYAqILdyd43ncPxT7A5O8bzuH4p9gCogt3J3jedw/FPsDk7xvO4fin2AKiC3cneN53D8U+wOTvG87h+KfYAqILdyd43ncPxT7A5O8bzuH4p9gCogt3J3jedw/FPsDk7xvO4fin2AL13SwFSrNODgloqL0r313urH0n/c+Xi6s1K0W03L3Jazsw8noy2talf8AC49q2KM1KDdNQve1vOd72d8zehTcY2dtew2i3a7VnbNa7HLgsRUqOcpRUYXtFbX09X+dLm6xFk84Tc0046OV073NpwvqM3lpW0fJt519u6xHWr2lGnHz5K66EtbM4sUYYqasYqg6kYrLJpu5JRhoxS3EWOqVI0pOkk55WTTa155LoI+5lavODdeChLSskk1dWWeftNU3SuzodPyrpvXmZlTu73djjoVKjxM4OqmoRhJx0Y/zad1vVrR9/STYqpUjKOirxeXmt2eu+TJ6wOhkdKnJPN3yRu7qO922bzFKTau1YoxCDUm28tn/ACJU/LUrm0JSbd42SeTve/T0HF3Yr1KdBypecurPJuzbyjd2V3+omKjtlFvUJK9jj7kYmpUoKVReVlqad/JTausnZtq63EXcjun9ZV0tl5JaoXvaDe2p9pbPmH0HF6KS1mKCkl5dr3eokI6M5SjeUHBv+VtNrrtlcDNOGiraTlrzlrNKcainno6NunSuSRbu7xslqd739mwjnKemkl5Od3b3EmaLDmxuG06kG6amknk1ck7n0XCMloqN5XSWraZ7o1Kkaa721FuUVpytaN3ra2rZ7TfCVnOGlLRvdryXeLs7XT3FHzX3Fn9f+td98n7Od9VtHqIO7PcCpiK1SpGcEp06UEne94VVNvVuPqYvGd7mk7Rja+k3rztoxW19HUdUJNxi2rNq7W7LUZw0iZo1j80+SYidoor/ANLvo7U7oKj3ucId70r6V8721WXQfYwOGdHD0qTs3Tpxi2tTsksvcdN8z52OrVI6Tpxc5XS0dLRyuru73LM0w6KkJN3jJx1fyt6r/wBzHe52ac27q3mvLpN8LJ+VnfUc/cXF4itSlLE0O8TU2lG97pan817LhUkaU02++Szd7KIjSmmv3krK2Wjr9pNiqrjDyb6TyjZXz32MYSu6lKM3GUW1qktF9dthURypNq2lLXe9nfU/frMOnU5yXCQd28biqPefq2H7/pTtPX5K35aut5Kx9DETlGE3COlJJuMd7tkiCKjFxTu3Jt3vZrYl+hBDD1FV03Ulo3b0fKtZqyX6/wBiTuXiK1WjGdan3ubbvGzW3J2eaFOvVeIlB07U0spZ56tup9RYvdn2i366oajJy4+vVpqLpQ023Z5N293zJ6spKEnGOlJJtRvrdslczE1mixiiZmG4OLuRia1Wip16Xep3a0bNZbHZ5ozHEVvrLpun+6Sup57t+rXlY3MUmiuwHB3XxVelTjKhR77JySazdlvsjtk3otpZ2va+22q5mt6Nes0ieWwOPuXiK1WnpVqfe5XatZq632Znv9X6xod7/d286z3b9WvKwx/4zSVnBMTMcOsHBiJyUnm0r9L/AARNhJSad9dllfbuMxiicXqy6QRKbcZX3ZOzWzczgdSV9b1pWs9VtdxixxhpUfUBpRfkrqNzSAAAAAAAAAAAAAAAAOJy6F7jaNRrVZew0AEnf5bx3+W8jAGuMxbhTlJ1FC38zSdj4vhifr1P4aOn6R+iT64/mRSiMeTHGClls8MT9ep/DQ8MT9ep/DRUwVyzo47WzwxP16n8NDwxP16n8NFTAM6OO1s8Mz9ep/DQ8MT9ep/DRUzIM6OO1r8MT9ep/DQ8MT9ep/DRVADOjjta/DM/Xqfw0PDM/Xqfw0VMAzo47WzwxP16n8NDwxP16n8NFTAM6OO1s8MT9ep/DQ8MT9ep/DRUwDOjjtbPDE/Xqfw0PDE/Xqfw0VMAzo47WzwxP16n8NDwxP16n8NFTAM6OO1s8MT9ep/DRq+6jf8A3tL4aKqAudHHa1x7rSWrG0l/40beGZ+vU/hoqQCZ0cdrZ4Yn69T+Gh4Zn69T+GipgXM6OO1t8Mz9ep/DRjwzP16n8NFTAM6OO1t8Mz9ep/DRjwxP16n8NFTAM6OO1t8MT9ep/DQ8Mz9ep/DRUgDOjjtbfDM/Xqfw0PDM/Xqfw0VIAzo47W3wzP16n8NDwzP16n8NFSAuZ0cdrb4Zn69T+Gh4Zn69T+GipALnRx2tT7qN68bS+GjMe60lqxtJf+NFUMgzo47WvwxP16n8NGnhR+uUvhoq4F0zo47Wtd2Jr/vqXw0PDM/Xqfw0VQC650cdrX4Zn69T+Gh4Zn69T+GiqHXh4Ue81HNSc1azWWjd5dbfyFzOjjtYPDM/Xqfw0PDM/Xqfw0fAxlGnGnTnGLjKWuOlpKNknr15qUXbZvOMXJ80Rt2tfhmfr1P4aHhmfr1P4aKoBczo47WvwxP16n8NH3aOJlKEXpKV1rSWfSebl67jei0vuojpgxRiitH0e/S3jv0t5GCtpO/S3jv0t5GAAAAAAD5n0j9En1x/MilF1+kfok+uP5kUoOHn2ZO7vUHGnFq0pwumltTnrz22ijjpaOktO+jtsdX8L/Wc8c/9ccLKwUftbUlq+xpf3XuMYSlDSqXWlFJWyz8+K1X3Nj+F/rH8L/Wc6zTf+LSCWDin598pvyc1aKllf/aveSQwcU5JtS2LVl5UPxs2R/wv9Y/hf6xXF+/xbNpYBaVlJ61s3qTy4fxI1hYupKOk0o21qzd2lez2Z3JISw0XdaafVv1mJSwz16YjFi/f4UhpVorSpwVrOGk5La7Nv5DvKi6iVpSjJKPSs7u23Z7zb+F/rH8L/WWMU/v8ShDCRnZ30dJxyVna83FrXsVmaxw0G4pTau3HO2TWv9LdZt/C7pj+F/rHtP7/AAow8PBU5tPS8nybpJpqcU112b9hK8NGVOGSjLRi76r+de2eb1Ef8L/WP4X+skzP7/BpisMqa85tuTSyVrJJ3/H8DlOqr9X0XoKWlsucp2wTMxdmQAG2QAAAZMAAZAGAZMAAAAAAAAAAAAAAAAAAAAADA663eYNR73J+TFt98tm4pvLR6TeM6PeZ/u5pacL2qLdPfEi7of6sramoNe2EWR0q2ipRcdKMrXV2tWrP2hqt3TWxNKpGKaqLR2pRbeSWerZFGtChRqS0IurpNPRvGNr2bV7NnS8PQpxp1Z03OnP7M5ZdLWimujPM5sFovEx0IvRcnoxvd2s7K+8LvdHHA120lSqZ/wBLRHWoypzcJrRlHJp7Drq01RgnKhOEpxvTbq+bb+aySd10keKwsaSSlJym76laKalaSu3d6nnbcEo5S9dxvRaX3UUUvXcb0Wl91E3ejwf6y7QAHYAAAAAAAB8z6R+iT64/mRSi6/SP0SfXH8yKUHDz7AAK8wAAAAAAAAAAAAAAAAAb06cpu0U5PciVorUwdDpQh58rv7MP1lqXsuY+tNeZFQ6VnLiefusZ9q6LTliOFm1fR0Vvk1Fe9me8wXnVV/ti5fi7IhlJt3bbe95swKTylk96K2VJf7ox/Cz+Zh1obKUfa5v5NEIL6wVT/WFzVP3S7Q+sLmqful2iEwT1gqmVaG2lD2Oa+bZm9F7KkeqUZfhZfMgBfWCqfvUH5tVf7ouP4q6MSws0rqOkt8WpL3rUQmYtp3Tae9ZMlJ5LMAn+tN+elPpfncSz99x3uEvNlov7M/0nq96QrTUogBvUpyg7STT6f03mprUYABUAAAAAAAAAAB0vF3tenSbsldqV3ZJLVLcjH1iO2jT9mmv/AGOcBaujv1KzXepJPdUa+aZvRr0YvSVOpfRkl+8jti1fzOk5ACojqp4+aTjK04uOjaS1LZZ7H0nKAVdOKxMamjalCFsvJul7i5dxvRaX3UUUvXcb0Wl91E3erw6S7QAHUAAAAAAAB8z6R+iT64/mRSi6/SP0SfXH8yKUHDz7AAK8wAAAAAAAAAAAAAGUjelSc3lqWbbySW9skdZQyp69s3lJ9X2V+PyMzO0LQ7zGH+pr+xF5/wC57PxfUYeKlla0Yp30Y5L2730shME9a6rWmiWtJytLf8yIkpv+V6n+D2P/ADeaNWLhtZcV/wDJgAGmAAAAAAAAAAADJgAS068orRylH7Ms1/w+lG/eoz/08pfYk/yvb1Oz6znBn13hastWyeTRg6FWUsqt3umvOXX9pfiR1aTg1ezTzi1qa6BE7SURgGTSMAAAAAAAAAAAAAAAAyXruN6LS+6iil67jei0vuom71eDSXaAA7AAAAAAAAPmfSP0SfXH8yKUXX6R+iT64/mRSg4efYABXmAAAAAAAyBg2pwcpJLayXvKz8r+W633tfV7yOnJp5bcveWlJuzWsWZrwUZOKd7GcPR03a9kk29rsty2sxWk23fezSMmmmnZrU0THvRrDpFUlWtdKMVowWpfq3tf+IjJ5xVROUUlJZzitv8AVH9V+mrnM4aLIADSBJPNKW3U+vYyMmowdruyg8m3kvZvfUZxWu3h4QmYxbdkm3uWZLLQg7W0mtssl7l/cxKvNq17Lcsl7kKzOiUpqd4lt0Y/ekk/drHe4bai9kZP52IjIpPKJNGn9qXAu0NGn9qfAu0RGRT9Ene4c574tfK47w9jjLqkr+55kZgUnkbTg4u0k0+lWMHTSraNNrflnmrrPU8t3vItOEvOjovfHsv9GjMYp4anDSl0ZgllRdrpqUVra2da1oiNxNWQAFQJaVbR8lrSg9cf1T2PpIgSYqqWrStaUXpQep/o1sZJSUdGzzck/ZbU/mZw3kxcp/6bycftvcurfs/Ago+fHrXzHjxUmjOOKw1krPXfqBJTpKSbbtb5bTWpFKTS1Jmpi1SMUVo0ABFAAAAAAAAAABkvXcb0Wl91FFL13G9FpfdRN3q8Gku0AB2AAAAAAAAfM+kfok+uP5kUouv0j9En1x/MilBw8+wACvMAAAAAAAA3U3dPdqE3HLRTWSv1mpgVKJK3nN78/eRki8qNtq/FGhZ5SOGYyaaadms00TTiqic4q0lnOK/NHo3rZ1auc2hNxaadms00YmN4aYMwg5OyV2TuEanlRcYP+dPJL+qPR0Ec6itowyjt3y6/7E9q6FGbxhunLf8Ayr+/y6yOc3J3k7s1BYgSSzjfasn1bH+nuIzaErPo1PqE42dhFpo1N4q1ABpgAAAA3pLO71LN+z/mxJmixFZozVyaj9lW9utmgbuYERSFxTWWYyad02mtTRLpRn51oy+0lk+tfqvcQgkxVG9Sm4uz9m59KZqbwqWWjJXju2rpW5ipTtZp3i9T/R7n0Cu0jQmp00lpz83+WO2X9lvYp01FKdRa/Nj9rpe6PzI6lRyd5PP5dCWxErW0LoVajm7vqSWSS2JLYjNPK73L8dSNErm88lor29ZvDFGJvZpfK2w3qOLSf8z1kZkVWjAAAAAAAAAAAAADJeu43otL7qKKXruN6LS+6ibvV4NJdoADsAAAAAAAA+Z9I/RJ9cfzIpRdfpH6JPrj+ZFKDh59gAFeYAAAAAAAAAAGU7ajeylqylu2Pq/sRgVJhlowSKpfKSv8/eNGP2vei04SvKMyZklsd/ZY1IoAAMm/nR6Y/L/j9TVRv7riErO5mWotqwYN6kbPLVrXUaFiapMUmgACoyb6ofef4L/n5GiJKkXfJNqOV9mWv8TM8N4YtMogAVgABQJKdVx2Jp60809xouk38j+r8BSpWjWc3JuUndvWzMYN56lveozppal7Xn/wayk3rdxSIS8ttJLKPtf9jQwBVYgAAAAAAAAAAAAAAABkvXcb0Wl91FFL13G9FpfdRN3q8Gku0AB2AAAAAAAAfM+kfok+uP5kUouv0j9En1x/MilBw8+wACvMAAAAAAAAAAAAAAAAAAAAAN6PnLpy9+X6m9WmorU1nbPb09BCZMzF6txitRvHyo22rNfqv195GZjJp3WtG1SOd1qea/sNJNYaAA0wkpZPS+yr+3Z+JuqsdGzT1JPotnl1s0eUUt+b9mS/UjMesYry6e04bQy3d6rGADTAACoHV9Rnkm4KUo6UYt+U1a62W1dJyn0qWIjorRdDzFCXfHNTSSzStsbvqz+QWHDRoynfRtZZttpJbrt5Gh1zdKUVCFRQjGcn5elmnbReSeaSe78TWGIp94dOd7p3ho21v7V9nUCjlAuZAwAAgAAAAAAAAAAAAAyXruN6LS+6iil67jei0vuom71eDSXaAA7AbtmDSr5kup/IDcAAAAB8z6R+iT64/mRSi6/SP0SfXH8yKUHDz7AAK8wAAAAAAAADJgAAAAAAAAAAAAAAybwzTj7V17vaRmUyTFWomksGYxbaS1vI3qK/lLU/we3/ADpFPK8tyy63/jJW1V9b0YqyvJ21LJdSNDJgsRRmZrNQAFQAMgE7E/16vztTiZ0wtCk5qELaCtOSUnpuSvFXyVlpZW1Z7T54XRP9er87U42Pr1fnanExTpwdNttaWdvLt/8AXRd/eT0tHvK/09HRm5X0dPTzUEv5l/K8stYW7mqYurJaMqk5J605NohB14arRjB6Sl3y7tJJPRyVnZuzzvl+KCauaUXF2as7J+xq6/BmLHd3WrOc4yUlKDjFxWV45JWlbbdPIgwWI71VjPPK+rXqaBS7RYapfR73PSte2i7232NY05OSil5Tdkum9rG3f2pymkvKbdpJTtd3/mT95rCo4yUk7STumtjA3r4dwno+dle6UtVr3zWrpISeWMqNWcteTdlpNXvZy12u9RACQABAAAAABkvXcb0Wl91FFL13G9FpfdRN3q8Gku0AB2DSr5kup/I3NKvmS6n8gNwAAAAHzPpH6JPrj+ZFKLr9I/RJ9cfzIpQcPPsAArzAAAAAAAb0qbnJRWt79S3t9FiTNFbQpXhKbdksl0yez3XZGS4ionaMfMjlHp3yfS/7ERMNdZJYABpAAAAAAAAAAAAABJSz8l7dXXs/sJ5RS9r/AE/zpIzaUm3du7M0u37Wo1ABpgAAAAyBvQnozi7tJSTdtdr5ndVxlNyrrRWhKNqSUIpJ/a3rb7yCnhI97c6k2slJRiruzbWbdktT3nKGrwwDanBydoq7ZPRw37xxqNxSi5ScbSaSV1ts9m3aEcwJK9LQdr3TSaeq6erLYa04aUlFO12ld6s2BqD7EMFRjpStpaNOd4uavpJO18tbTvZarPPI+Q9b2BZijAOyrKkqUJU4x05aSmpPScLaml0rO/uI61WEoJLSurbIJdOaV37QlESpvQctiaT9qbXyZodWDq00qkamqSVrptXTyvZp6rkWJcHNuCtGy2Wztm0ruyv0gRAAIAAAAAMl67jei0vuoopeu43otL7qJu9Xg0l2gAOwaVfMl1P5G5pV8yXU/kBuAAAAA+Z9I/RJ9cfzIpRdfpH6JPrj+ZFKDh59gAFeYAAAAACWFTRhJJZyyb223e39CIySYqoYMmCoAAAAAAAAAAAAAAAAAAAAAAAAGTBlJt2WbepIDpw9XS0KM3am5pNqyaTe/cm2yeNOjo15KC/dNKF5tqV5NXe/flb3GKvcyapxbSjLQk5LNt2cnsul5KWu1z54a01b056LvaL6JJSXuZvHESUtJKOas1opRaeyyyIjARJVqucru2pJJZJJKySIwANqc5QalFuMlqadmvaYMGQjAOmWCqKCm9GzjpJaS0tH7WjrscwVkGAEAAAAAAAAZL13G9FpfdRRS9dxvRaX3UTd6vBpLtAAdg0q+ZLqfyNzSr5kup/IDcAAAAB8z6R+iT64/mRSi6/SP0SfXH8yKUHDz7AAK8wAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAB0YN2c57YQk11u0U/fK5zk+H82t/wDGv/0gFgoxbjJuo4wjbS1vXklb2bTfFUacG4xlNyVtcUlqvvZBSqyg24vWrO6TTW5p5M+7UjQl5UoXloKSsruT0VeTSd5LSl0eaw1EVfAPpYbCtYWrOUVmouOpu2ee9b+pHN3QpwjUtDLJOSta0tqtd26rsgVSWXlPJNLN5J60gmjoxGEVOUINuU5Wfkp+a3aKSdnfXl7DeeHjB1YKUZrvd1JW1xkr23apI5ataU3ecnJ2td7iMFklHR0lpauvR/Gz+RviVBNaGja38spS97aX4IhMBHdOpC0pqSbdKMIxs7p6MYy2WtZS27ThABUAAQAAAAAAABkvXcb0Wl91FFL13G9FpfdRN3q8Gku0AB2DSr5kup/I3NKvmS6n8gNwAAAAHzPpH6JPrj+ZFKLr9I/RJ9cfzIpQcPPsAArzAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAHXKVJU3oN6Uoxi008ndSk76rXSSOU6qFSiqUozjLTk/OST0UtTXt1rLLaFhym8YNR74na0tHJ2d7X/wA6yTDwg6q1uCvLPW1FOVn12J6HdJxi04+U5aWlBxhsS0baLsstmsK4AbVJXk3ZRTbdlqXQjW4QAOmlQUqFSpfyoSgrdEtL+wHMCSFGUrWi3duKttazaXTmbRo+TO91KMVJLem0v/ZMCEGTAQAAAAAAAAAAGS9dxvRaX3UUUvXcb0Wl91E3erwaS7QAHYNKvmS6n8jc0q+ZLqfyA3AAAAAfM+kfok+uP5kUouv0j9En1x/MilBw8+wACvMAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAOjCf9R7qcvxtH/2M4PDKrLQ0nGbvbybxsld3d7317DanRlCnOUlaMoJRd1m3KLy9iZzQm4u8W01qadmGnXhEtKbpLSknHQUrOWjfynbfq1arvrN8VRr1HLRjKdNTkotRusm1sWZ8878L3TlShGKirxb0W7ZKTTlsvna2vbqBH65K9CdN6M4uL3NHbV7pJwcIRlGOuPlWcZLJWaWq1+l3buQY7F99cfOtFPzmm8229SX+I5QVpo7auNU50pSi5d7ilJOXnNNvS9uV+o1njNKdSco+VUg4ys9cm09O2zNajkASrJgAAAAgAAAAAAADJeu43otL7qKKXruN6LS+6ibvV4NJdoADsGlXzJdT+RuaVfMl1P5Abg6Pq63sfV1vYHODp+qdI+qdIHw/pH6JPrj+ZFKPTMT3MjVg4Tzi9Z81/RPCrZLiZGceCMdLqKC9eKmF3S4mPFTC7pe9lu55Ec9fVFBevFTC7pe9jxUwu6XvYuZEc9fVFBevFTC7pe9mPFXC7pcTFzIjnr6owLz4q4XdLiZnxUwu6XExcyI56+qKC9eKmF3S4mPFTC7pcTFzIjnr6ooL14qYXdLiY8VMLulxMXMiOevqigvXiphd0uJjxUwu6XExcyI56+qKC9eKmF3S4mPFTC7pcTFzIjnr6ooL14qYXdLiY8VMLulxMXMiOevqigvXiphd0uJjxUwu6XExcyI56+qKC9eKmF3S4mPFTC7pcTFzIjnr6ooL14qYXdLiY8VMLulxMXMiOevqigvXiphd0uJjxUwu6XExcyI56+qKC9eKmF3S4mPFTC7pcTFzIjnr6ooL14qYXdLiY8VMLulxMXMiOevqnYSVGOm6kXN2WjHUnrvd7Nm/aRVtC/7vStvla/uRdvFTC7pcTHiphd0uJi5kxz19UUF68VMLulxMeKmF3S4mLmRHPX1RQXrxUwu6XEx4qYXdLiYuZEc9fVFBevFTC7pcTHiphd0uJi5kRz19UUF68VMLulxMeKmF3S4mLmRHPX1RQXrxUwu6XEx4qYXdLiYuZEc9fVFBevFTC7pcTHinhd0uJi5kRz19UUF68U8LulxMeKmF3S97FzIjnr6ooL14qYXdL3seKmF3S4mLmRHPX1Ri9dxvRaX3UY8VMLulxM+nRwUYRUI5JKyI6YMEYIpVECaVGK1t+4wqcb20nfqK0iNKvmS6n8jr+rrezEsKmmm3mAqwqaV4SSySs9WvN/IYdVf+o10Ja9e1nnXKJjeaw/DPtjlExvNYfhn2yj02jGS0tJ65NrO9lsNpptZOz32ueYcouO5uhw1O2OUXHc3h+Gp2wPUGaVL7NdsjzLlFx3N0OGp2xyiY3m8Pw1O2B6LGnU2yu9JN56ktiyJJp6L0dew825RMbzWH4Z9scomN5rD8M+2B6E/rF3bQW7W9+XyFdV9L924W8nztlm9KyttVtuw895RMbzWH4Z9scomN5rD8M+2B6BhvrLknUdNRybST0lk7q+rXYloqSilJ6Ulre8865RMbzWH4Z9s1f7QcY/8ApYfhqdsD0pSdlfJ5XsSHmK/aDjOaw/DU7ZtyiY3msPwz7YHpgPM+UTG81h+GfbHKJjeaw/DPtgemA8z5RMbzWH4Z9scomN5rD8M+2B6YDzPlExvNYfhn2xyiY3msPwz7YHpgPM+UTG81h+GfbHKJjeaw/DPtgemA8z5RMbzWH4Z9scomN5rD8M+2B6YDzPlExvNYfhn2xyiY3msPwz7YHpgPM+UTG81h+GfbHKJjeaw/DPtgemA8z5RMbzWH4Z9scomN5rD8M+2B6YDzPlExvNYfhn2xyiY3msPwz7YHpgPM+UTG81h+GfbHKJjeaw/DPtgemA8z5RMbzWH4Z9scomN5rD8M+2B6YDzPlExvNYfhn2xyiY3msPwz7YHpgPM+UTG81h+GfbHKJjeaw/DPtgemA8z5RMbzWH4Z9scomN5rD8M+2B6YDzPlExvNYfhn2xyiY3msPwz7YHpgPM+UTG81h+GfbHKJjeaw/DPtgemA8z5RMbzWH4Z9scomN5rD8M+2B6RNSvkR4ylKVFxik5Xi7N2TtJNxv0pNHnfKJjeaw/DPtjlExvNYfhn2wPQe51OaU3KmqWlJNQUlK3kpa1lrJkpaXR7DzjlExvNYfhn2xyiY3msPwz7YHo2IjJ20em+diR30emx5ryiY3msPwz7Y5RMbzWH4Z9sD0mlpWel+n6EVTFxi2tCq7bqcmvfY875RMbzWH4Z9scomN5rD8M+2IHoccVGTS0KivvpyS9rtkTSSvqvbPJHmr/aHjeaw/DPthftDxq/6WH4anbEj0x6siOipK+l+n6HnHKJjeaw/DPtjlExvNYfhn2wKiAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAD//Z\n",
      "text/html": [
       "\n",
       "        <iframe\n",
       "            width=\"400\"\n",
       "            height=\"300\"\n",
       "            src=\"https://www.youtube.com/embed/tp71tZOeitk\"\n",
       "            frameborder=\"0\"\n",
       "            allowfullscreen\n",
       "        ></iframe>\n",
       "        "
      ],
      "text/plain": [
       "<IPython.lib.display.YouTubeVideo at 0x1e502c987c0>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.display import YouTubeVideo\n",
    "\n",
    "YouTubeVideo('tp71tZOeitk')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this file, we are going to consider the problem of the pure bending of a curved bar. And, again we are going to consider a $\\theta$-independent Airy stress function. \n",
    "\n",
    "We begin by importing the `polarUtilities`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from polarUtilities import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now set the Airy stress function to be $\\theta$-independent:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\phi{\\left(r \\right)}$"
      ],
      "text/plain": [
       "phi(r)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "phi = sym.Function('phi')(r)\n",
    "phi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now solve the biharmonic equation $\\nabla^4 \\phi = 0$ as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "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": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sym.dsolve(polarbiharmonic(phi))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This solution is, of course, the same solution that we had obtained in the first file. However, in order to have the same notation as in Timoshenko and Goodier, we rewrite $\\phi$ as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle A \\log{\\left(r \\right)} + B r^{2} \\log{\\left(r \\right)} + C r^{2} + D$"
      ],
      "text/plain": [
       "A*log(r) + B*r**2*log(r) + C*r**2 + D"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "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",
    "phi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We again use a better way to display the variables as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import Math, Latex"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\phi = 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": [
    "display(Math(r'\\phi = {}'.format(sym.latex(phi))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now obtain the expressions of the stress components $\\sigma_{rr}$, $\\sigma_{\\theta\\theta}$, and $\\sigma_{r\\theta}$:"
   ]
  },
  {
   "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": [
    "We note that up to this point there is no difference with the first file. However, we need to solve for $A$, $B$, $C$, and $D$ subject to the boundary conditions appropriate for this problem, which are:\n",
    "\n",
    "\\begin{align}\n",
    "\\sigma_{rr} &= 0 \\quad \\text{at $r=a$} \\\\\n",
    "\\sigma_{rr} &= 0 \\quad \\text{at $r=b$} \\\\\n",
    "\\int_a^b \\sigma_{\\theta\\theta} \\; {\\rm d}r &= 0 \\\\\n",
    "\\int_a^b r \\sigma_{\\theta\\theta} \\; {\\rm d}r &= -M_0\n",
    "\\end{align}\n",
    "where $M_0$ is the externally applied moment. \n",
    "\n",
    "We now set up these boundary conditions as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{A}{a^{2}} + 2 B \\log{\\left(a \\right)} + B + 2 C = 0$"
      ],
      "text/plain": [
       "Eq(A/a**2 + 2*B*log(a) + B + 2*C, 0)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{A}{b^{2}} + 2 B \\log{\\left(b \\right)} + B + 2 C = 0$"
      ],
      "text/plain": [
       "Eq(A/b**2 + 2*B*log(b) + B + 2*C, 0)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{A}{b} - \\frac{A}{a} - 2 B a \\log{\\left(a \\right)} + 2 B b \\log{\\left(b \\right)} - a \\left(B + 2 C\\right) + b \\left(B + 2 C\\right) = 0$"
      ],
      "text/plain": [
       "Eq(A/b - A/a - 2*B*a*log(a) + 2*B*b*log(b) - a*(B + 2*C) + b*(B + 2*C), 0)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle A \\log{\\left(a \\right)} - A \\log{\\left(b \\right)} - B a^{2} \\log{\\left(a \\right)} + B b^{2} \\log{\\left(b \\right)} - a^{2} \\left(B + C\\right) + b^{2} \\left(B + C\\right) = - M_{0}$"
      ],
      "text/plain": [
       "Eq(A*log(a) - A*log(b) - B*a**2*log(a) + B*b**2*log(b) - a**2*(B + C) + b**2*(B + C), -M_0)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "a, b = sym.symbols('a, b', positive=True)\n",
    "M0 = sym.symbols('M_0')\n",
    "\n",
    "lhs1 = sigmarr.subs(r,a)\n",
    "rhs1 = 0\n",
    "eq1 = sym.Eq(lhs1,rhs1)\n",
    "display(eq1)\n",
    "\n",
    "lhs2 = sigmarr.subs(r,b)\n",
    "rhs2 = 0\n",
    "eq2 = sym.Eq(lhs2, rhs2)\n",
    "display(eq2)\n",
    "\n",
    "lhs3 = sym.integrate(sigmatt,(r,a,b))\n",
    "rhs3 = 0\n",
    "eq3 = sym.Eq(lhs3,rhs3)\n",
    "display(eq3)\n",
    "\n",
    "lhs4 = sym.integrate(r*sigmatt,(r,a,b))\n",
    "rhs4 = -M0\n",
    "eq4 = sym.Eq(lhs4,rhs4)\n",
    "display(eq4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we solve for $A$, $B$, and $C$. Note that while we have only three unknowns ($D$ does not appear), we have four equations. Strictly, speaking, we should be using only three of these equations to solve for the three unknowns (one of the equations is not independent). However, operationally speaking, we don't need to worry. A naive implementation using `linsolve` using all four equations still gives us the correct solution; thus:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle A = \\frac{M_{0} \\log{\\left(\\left(\\frac{a}{b}\\right)^{4 a^{2} b^{2}} \\right)}}{a^{4} - 4 a^{2} b^{2} \\log{\\left(a \\right)}^{2} - 4 a^{2} b^{2} \\log{\\left(b \\right)}^{2} - 2 a^{2} b^{2} + b^{4} + \\log{\\left(b^{\\log{\\left(a^{8 a^{2} b^{2}} \\right)}} \\right)}}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle B = \\frac{2 M_{0} \\left(a^{2} - b^{2}\\right)}{a^{4} - 4 a^{2} b^{2} \\log{\\left(a \\right)}^{2} + 8 a^{2} b^{2} \\log{\\left(a \\right)} \\log{\\left(b \\right)} - 4 a^{2} b^{2} \\log{\\left(b \\right)}^{2} - 2 a^{2} b^{2} + b^{4}}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle C = \\frac{M_{0} \\left(- a^{2} + b^{2} + \\log{\\left(a^{- 2 a^{2}} b^{2 b^{2}} \\right)}\\right)}{a^{4} - 4 a^{2} b^{2} \\log{\\left(a \\right)}^{2} - 4 a^{2} b^{2} \\log{\\left(b \\right)}^{2} - 2 a^{2} b^{2} + b^{4} + \\log{\\left(b^{\\log{\\left(a^{8 a^{2} b^{2}} \\right)}} \\right)}}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "soln, = sym.linsolve([eq1,eq2,eq3,eq4],(A,B,C))\n",
    "\n",
    "Asoln = soln[0]\n",
    "Bsoln = soln[1]\n",
    "Csoln = soln[2]\n",
    "\n",
    "display(Math(r'A = {}'.format(sym.latex(Asoln))))\n",
    "display(Math(r'B = {}'.format(sym.latex(Bsoln))))\n",
    "display(Math(r'C = {}'.format(sym.latex(Csoln))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have the values of $A$, $B$, and $C$, we can make the appropriate substitutions to find the solutions to the stress components as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\sigma_{rr}^{\\rm soln} = - \\frac{4 M_{0} \\left(- a^{2} r^{2} \\log{\\left(r \\right)} + b^{2} r^{2} \\log{\\left(r \\right)} + r^{2} \\log{\\left(a^{a^{2}} b^{- b^{2}} \\right)} + \\log{\\left(a^{- a^{2} b^{2}} b^{a^{2} b^{2}} \\right)}\\right)}{r^{2} \\left(a^{4} - 4 a^{2} b^{2} \\log{\\left(a \\right)}^{2} - 4 a^{2} b^{2} \\log{\\left(b \\right)}^{2} - 2 a^{2} b^{2} + b^{4} + \\log{\\left(b^{\\log{\\left(a^{8 a^{2} b^{2}} \\right)}} \\right)}\\right)}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\sigma_{\\theta\\theta}^{\\rm soln} = - \\frac{4 M_{0} \\left(- a^{2} r^{2} \\log{\\left(r \\right)} - a^{2} r^{2} + b^{2} r^{2} \\log{\\left(r \\right)} + b^{2} r^{2} + r^{2} \\log{\\left(a^{a^{2}} b^{- b^{2}} \\right)} + \\log{\\left(a^{a^{2} b^{2}} b^{- a^{2} b^{2}} \\right)}\\right)}{r^{2} \\left(a^{4} - 4 a^{2} b^{2} \\log{\\left(a \\right)}^{2} - 4 a^{2} b^{2} \\log{\\left(b \\right)}^{2} - 2 a^{2} b^{2} + b^{4} + \\log{\\left(b^{\\log{\\left(a^{8 a^{2} b^{2}} \\right)}} \\right)}\\right)}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sigmarr_soln = sigmarr.subs([(A,Asoln),(B,Bsoln),(C,Csoln)]).simplify()\n",
    "sigmatt_soln = sigmatt.subs([(A,Asoln),(B,Bsoln),(C,Csoln)]).simplify()\n",
    "\n",
    "display(Math(r'\\sigma_{{rr}}^{{\\rm soln}} = {}'.format(sym.latex(sigmarr_soln))))\n",
    "display(Math(r'\\sigma_{{\\theta\\theta}}^{{\\rm soln}} = {}'.format(sym.latex(sigmatt_soln))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note even though these solutions look quite different from those given in Timoshenko and Goodier, they are exactly the same. They may be easily simplified through elementary algebraic identities. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we will find out the displacements corresponding to the stress fields. For this, we will first find out the strains using the generalized Hooke's law (linear, elastic, isotropic material behaviour), and then find the displacements using the strain-displacement relations. \n",
    "\n",
    "For the application of the Hooke's law, we use the `sigmarr` and `sigmatt` variables instead of the `sigmarr_soln` and `sigmatt_soln` variables because the expressions of the former pair of variables contain $A$, $B$, and $C$ without their actual values substituted. Thus, we have:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "E, nu = sym.symbols('E, nu')\n",
    "\n",
    "ep_rr_stress = 1/E*(sigmarr - nu*sigmatt)\n",
    "ep_tt_stress = 1/E*(sigmatt - nu*sigmarr)\n",
    "ep_rt_stress = (1+nu)/E*sigmart"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\varepsilon_{r\\theta} = \\frac{\\frac{A}{r^{2}} + 2 B \\log{\\left(r \\right)} + B + 2 C - \\nu \\left(- \\frac{A}{r^{2}} + 2 B \\log{\\left(r \\right)} + 3 B + 2 C\\right)}{E}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\varepsilon_{\\theta\\theta} = \\frac{- \\frac{A}{r^{2}} + 2 B \\log{\\left(r \\right)} + 3 B + 2 C - \\nu \\left(\\frac{A}{r^{2}} + 2 B \\log{\\left(r \\right)} + B + 2 C\\right)}{E}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\varepsilon_{r\\theta} = 0$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display(Math(r'\\varepsilon_{{r\\theta}} = {}'.format(sym.latex(ep_rr_stress))))\n",
    "display(Math(r'\\varepsilon_{{\\theta\\theta}} = {}'.format(sym.latex(ep_tt_stress))))\n",
    "display(Math(r'\\varepsilon_{{r\\theta}} = {}'.format(sym.latex(ep_rt_stress))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The strain-displacement relations in polar coordinates are:\n",
    "\n",
    "\\begin{align*}\n",
    "\\varepsilon_{rr} &= \\frac{\\partial u}{\\partial r}, \\\\\n",
    "\\varepsilon_{\\theta\\theta} &= \\frac{1}{r} \\frac{\\partial v}{\\partial \\theta} + \\frac{u}{r}, \\\\\n",
    "\\varepsilon_{r\\theta} &= \\frac{1}{2} \\left( \\frac{\\partial v}{\\partial r} - \\frac{v}{r} + \\frac{1}{r}\\frac{\\partial u}{\\partial \\theta} \\right),\n",
    "\\end{align*}\n",
    "where $u(r, \\theta)$ and $v(r,\\theta)$ are the displacement field variables along the radial and the azimuthal directions, respectively. Note that each of them is a function of both $r$ and $\\theta$.\n",
    "\n",
    "We first consider the first of the strain-displacement relations involving $\\varepsilon_{rr}$, to determine the form of $u(r,\\theta)$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle u{\\left(r,\\theta \\right)} = \\frac{- A \\nu - A - 2 B r^{2} \\left(\\nu - 1\\right) \\log{\\left(r \\right)} + E r \\operatorname{f_{1}}{\\left(- \\theta \\right)} + r^{2} \\left(- B \\nu - B - 2 C \\nu + 2 C\\right)}{E r}$"
      ],
      "text/plain": [
       "Eq(u(r, theta), (-A*nu - A - 2*B*r**2*(nu - 1)*log(r) + E*r*f_1(-theta) + r**2*(-B*nu - B - 2*C*nu + 2*C))/(E*r))"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "u = sym.Function('u')(r,theta)\n",
    "lhs = sym.diff(u,r)\n",
    "rhs = ep_rr_stress\n",
    "eq1 = sym.Eq(lhs,rhs)\n",
    "f1 = sym.Function('f_1')\n",
    "sol = sym.pdsolve(eq1,solvefun=f1)\n",
    "sol"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We rewrite the expression on the right hand side of this solution in a slightly more convenient form and assign it to the variable `u_interim`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle A \\left(- \\frac{\\nu}{E r} - \\frac{1}{E r}\\right) + B \\left(- \\frac{2 \\nu r \\log{\\left(r \\right)}}{E} - \\frac{\\nu r}{E} + \\frac{2 r \\log{\\left(r \\right)}}{E} - \\frac{r}{E}\\right) + C \\left(- \\frac{2 \\nu r}{E} + \\frac{2 r}{E}\\right) + \\operatorname{f_{1}}{\\left(- \\theta \\right)}$"
      ],
      "text/plain": [
       "A*(-nu/(E*r) - 1/(E*r)) + B*(-2*nu*r*log(r)/E - nu*r/E + 2*r*log(r)/E - r/E) + C*(-2*nu*r/E + 2*r/E) + f_1(-theta)"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "u_interim = sol.rhs.expand().collect(A).collect(B).collect(C)\n",
    "u_interim "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We next consider the second of the strain-displacement relations involving $\\varepsilon_{\\theta\\theta}$ to determine the form of $v(r,\\theta)$. Note that we make use of the just-determined form of $u(r,\\theta)$ (stored in the variable `u_interim`) to determine $v(r,theta)$. We store the solution in the variable `v_interim`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{4 B r \\theta}{E} + \\operatorname{f_{2}}{\\left(r \\right)} - \\int\\limits^{\\theta} \\operatorname{f_{1}}{\\left(- \\xi \\right)}\\, d\\xi$"
      ],
      "text/plain": [
       "4*B*r*theta/E + f_2(r) - Integral(f_1(-xi), (xi, theta))"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v = sym.Function('v')(r,theta)\n",
    "lhs = sym.diff(v,theta)\n",
    "rhs = r*ep_tt_stress - u_interim\n",
    "f2 = sym.Function('f_2')\n",
    "eq2 = sym.Eq(lhs,rhs)\n",
    "sol = sym.pdsolve(eq2,solvefun=f2)\n",
    "v_interim = sol.rhs\n",
    "v_interim"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we substitute the forms of $u(r,\\theta)$ and $v(r,\\theta)$ (stored, respectively, in `u_interim` and `v_interim`) in the final strain-displacement relation involving $\\varepsilon_{r\\theta}$ to obtain an equation involving both the unknown functions $f_1(\\theta)$ and $f_2(r)$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{4 B \\theta}{E} + \\frac{d}{d r} \\operatorname{f_{2}}{\\left(r \\right)} - \\frac{\\frac{4 B r \\theta}{E} + \\operatorname{f_{2}}{\\left(r \\right)} - \\int\\limits^{\\theta} \\operatorname{f_{1}}{\\left(- \\xi \\right)}\\, d\\xi}{r} - \\frac{\\left. \\frac{d}{d \\xi_{1}} \\operatorname{f_{1}}{\\left(\\xi_{1} \\right)} \\right|_{\\substack{ \\xi_{1}=- \\theta }}}{r} = 0$"
      ],
      "text/plain": [
       "Eq(4*B*theta/E + Derivative(f_2(r), r) - (4*B*r*theta/E + f_2(r) - Integral(f_1(-xi), (xi, theta)))/r - Subs(Derivative(f_1(_xi_1), _xi_1), _xi_1, -theta)/r, 0)"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lhs = (sym.diff(v_interim,r) - v_interim/r + 1/r*sym.diff(u_interim,theta))\n",
    "rhs = 0\n",
    "eq = sym.Eq(lhs,rhs)\n",
    "eq"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This equation can be easily simplified:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{r \\frac{d}{d r} \\operatorname{f_{2}}{\\left(r \\right)} - \\operatorname{f_{2}}{\\left(r \\right)} + \\int\\limits^{\\theta} \\operatorname{f_{1}}{\\left(- \\xi \\right)}\\, d\\xi - \\left. \\frac{d}{d \\xi_{1}} \\operatorname{f_{1}}{\\left(\\xi_{1} \\right)} \\right|_{\\substack{ \\xi_{1}=- \\theta }}}{r} = 0$"
      ],
      "text/plain": [
       "Eq((r*Derivative(f_2(r), r) - f_2(r) + Integral(f_1(-xi), (xi, theta)) - Subs(Derivative(f_1(_xi_1), _xi_1), _xi_1, -theta))/r, 0)"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eq.simplify()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From this equation we extract the following two equations:\n",
    "\n",
    "\\begin{align*}\n",
    "r\\frac{d {\\rm f}_2(r)}{d r} - {\\rm f}_2(r) &= C_0, \\\\\n",
    "\\int {\\rm f}_1(-\\theta) \\; d \\theta + \\frac{d {\\rm f}_1(-\\theta)}{d \\theta} &= -C_0, \n",
    "\\end{align*}\n",
    "where $C_0$ is some constant. The presence of this constant makes the equations slightly different from those present in Timoshenko and Goodier. \n",
    "\n",
    "Note that in order to arrive at the form of the second equation from the previous cell, we have replaced $\\xi$ with $\\theta$ and $\\xi_1$ with $-\\theta$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider the first equation. We replace ${\\rm f}_2(r)$ by ${\\rm f}_{\\rm r}(r)$ so that this first equation becomes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle r \\frac{d}{d r} \\operatorname{f_{r}}{\\left(r \\right)} - \\operatorname{f_{r}}{\\left(r \\right)} = C_{0}$"
      ],
      "text/plain": [
       "Eq(r*Derivative(f_r(r), r) - f_r(r), C_0)"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "C0 = sym.symbols('C_0')\n",
    "fr = sym.Function('f_r')(r)\n",
    "lhs = r*sym.diff(fr,r) - fr\n",
    "rhs = C0\n",
    "eq = sym.Eq(lhs,rhs)\n",
    "eq"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " The solution to this equation is given by: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\operatorname{f_{r}}{\\left(r \\right)} = C_{1} r - C_{0}$"
      ],
      "text/plain": [
       "Eq(f_r(r), C1*r - C_0)"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sol = sym.dsolve(eq)\n",
    "sol"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use $F$ instead of $C_1$; thus:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - C_{0} + F r$"
      ],
      "text/plain": [
       "-C_0 + F*r"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F = sym.symbols('F')\n",
    "fr = -C0 + F*r\n",
    "fr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we consider the second equation. We replace ${\\rm f}_1(-\\theta)$ by ${\\rm f}_\\theta(\\theta)$. This is an integro-differential equation and in order to solve it, we first differentiate it with respect to $\\theta$ to obtain the following form:\n",
    "\n",
    "\\begin{gather*}\n",
    "{\\rm f}_\\theta + \\frac{d^2 {\\rm f}_\\theta}{d \\theta^2} = 0\n",
    "\\end{gather*}\n",
    "\n",
    "Note that as a result of this differentiation, the constant $C_0$ vanishes from this second equation. We set up this equation as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\operatorname{f_{\\theta}}{\\left(\\theta \\right)} + \\frac{d^{2}}{d \\theta^{2}} \\operatorname{f_{\\theta}}{\\left(\\theta \\right)} = 0$"
      ],
      "text/plain": [
       "Eq(f_\\theta(theta) + Derivative(f_\\theta(theta), (theta, 2)), 0)"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ft = sym.Function(r'f_\\theta')(theta)\n",
    "lhs = ft + sym.diff(ft,theta,2)\n",
    "rhs = 0\n",
    "eq = sym.Eq(lhs,rhs)\n",
    "eq"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The solution to this equation is given by:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\operatorname{f_{\\theta}}{\\left(\\theta \\right)} = C_{1} \\sin{\\left(\\theta \\right)} + C_{2} \\cos{\\left(\\theta \\right)}$"
      ],
      "text/plain": [
       "Eq(f_\\theta(theta), C1*sin(theta) + C2*cos(theta))"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sol = sym.dsolve(eq)\n",
    "sol"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, we use $H$ and $K$ instead of $C_1$ and $C_2$, respectively to maintain the same notation as in Timoshenko and Gooider; thus:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle H \\sin{\\left(\\theta \\right)} + K \\cos{\\left(\\theta \\right)}$"
      ],
      "text/plain": [
       "H*sin(theta) + K*cos(theta)"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "H, K = sym.symbols('H, K')\n",
    "ft = H*sym.sin(theta) + K*sym.cos(theta)\n",
    "ft"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The expression of $u(r,\\theta)$ then becomes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle A \\left(- \\frac{\\nu}{E r} - \\frac{1}{E r}\\right) + B \\left(- \\frac{2 \\nu r \\log{\\left(r \\right)}}{E} - \\frac{\\nu r}{E} + \\frac{2 r \\log{\\left(r \\right)}}{E} - \\frac{r}{E}\\right) + C \\left(- \\frac{2 \\nu r}{E} + \\frac{2 r}{E}\\right) + H \\sin{\\left(\\theta \\right)} + K \\cos{\\left(\\theta \\right)}$"
      ],
      "text/plain": [
       "A*(-nu/(E*r) - 1/(E*r)) + B*(-2*nu*r*log(r)/E - nu*r/E + 2*r*log(r)/E - r/E) + C*(-2*nu*r/E + 2*r/E) + H*sin(theta) + K*cos(theta)"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "u_semi = u_interim.subs(f1(-theta),ft)\n",
    "u_semi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And, the expression of $v(r,\\theta)$ becomes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{4 B r \\theta}{E} - C_{0} + F r + H \\cos{\\left(\\theta \\right)} - K \\sin{\\left(\\theta \\right)}$"
      ],
      "text/plain": [
       "4*B*r*theta/E - C_0 + F*r + H*cos(theta) - K*sin(theta)"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v_semi = 4*B*r*theta/E + fr - sym.integrate(ft,theta)\n",
    "v_semi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the solution form for $v(r,\\theta)$ given by `v_semi` is different from the one presented by Timoshenko and Goodier because of the presence of the $C_0$ term. This constant represents a rigid body translation term and if we neglect it then we can proceed exactly like Timoshenko and Goodier to determine the values of the rest of the constants, ie. the values of $F$, $H$, and $K$ by using the following displacement boundary conditions:\n",
    "\n",
    "\\begin{align*}\n",
    "u &= 0 \\quad \\text{at $r=r_0$}, \\\\ \n",
    "v &= 0 \\quad \\text{at $r=r_0$}, \\\\\n",
    "\\frac{\\partial v}{\\partial r} &= 0 \\quad \\text{at $r=r_0$ and $\\theta=0$},\n",
    "\\end{align*}\n",
    "where $\\displaystyle r_0 = \\frac{a+b}{2}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle A \\left(- \\frac{\\nu}{E r_{0}} - \\frac{1}{E r_{0}}\\right) + B \\left(- \\frac{2 \\nu r_{0} \\log{\\left(r_{0} \\right)}}{E} - \\frac{\\nu r_{0}}{E} + \\frac{2 r_{0} \\log{\\left(r_{0} \\right)}}{E} - \\frac{r_{0}}{E}\\right) + C \\left(- \\frac{2 \\nu r_{0}}{E} + \\frac{2 r_{0}}{E}\\right) + K = 0$"
      ],
      "text/plain": [
       "Eq(A*(-nu/(E*r_0) - 1/(E*r_0)) + B*(-2*nu*r_0*log(r_0)/E - nu*r_0/E + 2*r_0*log(r_0)/E - r_0/E) + C*(-2*nu*r_0/E + 2*r_0/E) + K, 0)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle F r_{0} + H = 0$"
      ],
      "text/plain": [
       "Eq(F*r_0 + H, 0)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle F = 0$"
      ],
      "text/plain": [
       "Eq(F, 0)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "r0 = sym.symbols('r_0')\n",
    "\n",
    "lhs = u_semi.subs([(theta,0),(r,r0)])\n",
    "rhs = 0\n",
    "eq1 = sym.Eq(lhs,rhs)\n",
    "\n",
    "lhs = v_semi.subs([(C0,0),(theta,0),(r,r0)])\n",
    "rhs = 0\n",
    "eq2 = sym.Eq(lhs,rhs)\n",
    "\n",
    "\n",
    "lhs = sym.diff(v_semi,r).subs([(theta,0),(r,r0)])\n",
    "rhs = 0\n",
    "eq3 = sym.Eq(lhs,rhs)\n",
    "\n",
    "display(eq1,eq2,eq3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 0$"
      ],
      "text/plain": [
       "0"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 0$"
      ],
      "text/plain": [
       "0"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{A \\nu + A + 2 B \\nu r_{0}^{2} \\log{\\left(r_{0} \\right)} + B \\nu r_{0}^{2} - 2 B r_{0}^{2} \\log{\\left(r_{0} \\right)} + B r_{0}^{2} + 2 C \\nu r_{0}^{2} - 2 C r_{0}^{2}}{E r_{0}}$"
      ],
      "text/plain": [
       "(A*nu + A + 2*B*nu*r_0**2*log(r_0) + B*nu*r_0**2 - 2*B*r_0**2*log(r_0) + B*r_0**2 + 2*C*nu*r_0**2 - 2*C*r_0**2)/(E*r_0)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "soln, = sym.linsolve([eq1,eq2,eq3],(F,H,K))\n",
    "\n",
    "Fsoln = soln[0]\n",
    "Hsoln = soln[1]\n",
    "Ksoln = soln[2]\n",
    "display(Fsoln,Hsoln,Ksoln)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the above solution, $\\displaystyle r_0 = \\frac{a+b}{2}$. Since the solutions of $A$, $B$, and $C$ were found earlier, `Fsoln`, `Hsoln`, and `Ksoln` are actually expressed in terms of known values. "
   ]
  }
 ],
 "metadata": {
  "hide_input": false,
  "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"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}