{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Strain and stress tensors in spherical coordinates\n", "\n", "This worksheet demonstrates a few capabilities of [SageManifolds](http://sagemanifolds.obspm.fr/) (version 1.0, as included in SageMath 7.5) in computations regarding elasticity theory in Cartesian coordinates.\n", "\n", "Click [here](https://raw.githubusercontent.com/sagemanifolds/SageManifolds/master/Worksheets/v1.0/SM_elasticity_spher.ipynb) to download the worksheet file (ipynb format). To run it, you must start SageMath with the Jupyter notebook, via the command sage -n jupyter" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*NB:* a version of SageMath at least equal to 7.5 is required to run this worksheet:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "'SageMath version 7.5.1, Release Date: 2017-01-15'" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "version()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First we set up the notebook to display mathematical objects using LaTeX rendering:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true }, "outputs": [], "source": [ "%display latex" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Euclidean 3-space and spherical coordinates" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We introduce the Euclidean space as a 3-dimensional differentiable manifold:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "3-dimensional differentiable manifold M\n" ] } ], "source": [ "M = Manifold(3, 'M', start_index=1)\n", "print(M)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We shall make use of spherical coordinates $(r,\\theta,\\phi)$:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Chart (M, (r, th, ph))\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "Chart (M, (r, th, ph))" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "spher. = M.chart(r'r:(0,+oo) th:(0,pi):\\theta ph:(0,2*pi):\\phi')\n", "print(spher)\n", "spher" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Spherical coordinates do not form a regular coordinate system of the Euclidean space. So declaring that they span $M$ means that, strictly speaking, the manifold $M$ is not the whole Euclidean space, but the Euclidean space minus some half plane (the azimuthal origin). However, in this worksheet, this difference will not matter. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The natural vector frame of spherical coordinates is" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "Coordinate frame (M, (d/dr,d/dth,d/dph))" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "spher.frame()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We shall expand vector and tensor fields on the orthonormal frame $(e_1, e_2, e_3)$ associated with spherical coordinates, which is related to the natural frame $(\\partial/\\partial r, \\partial/\\partial\\theta, \\partial/\\partial\\phi)$ displayed above by means of the following field of automorphisms: " ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "d/dr*dr + 1/r d/dth*dth + 1/(r*sin(th)) d/dph*dph" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "to_orthonormal = M.automorphism_field()\n", "to_orthonormal[1,1] = 1\n", "to_orthonormal[2,2] = 1/r\n", "to_orthonormal[3,3] = 1/(r*sin(th))\n", "to_orthonormal.display()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In other words, the change-of-basis matrix is " ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "[ 1 0 0]\n", "[ 0 1/r 0]\n", "[ 0 0 1/(r*sin(th))]" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "to_orthonormal[:]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We construct the orthonormal frame from the natural frame of spherical coordinates by this change of basis: " ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "Vector frame (M, (e_1,e_2,e_3))" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "e = spher.frame().new_frame(to_orthonormal, 'e')\n", "e" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "e_1 = d/dr" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "e[1].display()" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "e_2 = 1/r d/dth" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "e[2].display()" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "e_3 = 1/(r*sin(th)) d/dph" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "e[3].display()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "At this stage, the default vector frame on $M$ is the first one introduced, namely the natural frame of spherical coordinates: " ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "Coordinate frame (M, (d/dr,d/dth,d/dph))" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "M.default_frame()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Since we prefer the orthonormal frame, we declare" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false }, "outputs": [], "source": [ "M.set_default_frame(e)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then, by default, all vector and tensor fields are displayed with respect to that frame:" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "e_3 = e_3" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "e[3].display()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To get the same output as in Out[10], one should specify the frame for display, since this is no longer the default one:" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "e_3 = 1/(r*sin(th)) d/dph" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "e[3].display(spher.frame())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Displacement vector and strain tensor\n", "\n", "Let us define the **displacement vector** $U$ in terms of its components w.r.t. the orthonormal spherical frame:" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "U = U_1(r, th, ph) e_1 + U_2(r, th, ph) e_2 + U_3(r, th, ph) e_3" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "U = M.vector_field(name='U')\n", "U[:] = [function('U_1')(r,th,ph), function('U_2')(r,th,ph),\n", " function('U_3')(r,th,ph)]\n", "U.display()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The following computations will involve the metric $g$ of the Euclidean space. At the current stage of SageManifolds, we need to introduce it explicitly, as a Riemannian metric on the manifold $M$ (in a future version of SageManifolds, one shall to declare $M$ as an Euclidean space, and not merely as a manifold, so that it will come equipped with $g$):" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Riemannian metric g on the 3-dimensional differentiable manifold M\n" ] } ], "source": [ "g = M.riemannian_metric('g')\n", "print(g)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Since $e$ is supposed to be an orthonormal frame, we declare that the components of $g$ with respect to it are\n", "$\\mathrm{diag}(1,1,1)$:" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "g = e^1*e^1 + e^2*e^2 + e^3*e^3" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "g[1,1], g[2,2], g[3,3] = 1, 1, 1\n", "g.display()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The expression of $g$ with respect to the natural frame of spherical coordinates is then " ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "g = dr*dr + r^2 dth*dth + r^2*sin(th)^2 dph*dph" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "g.display(spher.frame())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The covariant derivative operator $\\nabla$ is introduced as the (Levi-Civita) connection associated with $g$: " ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Levi-Civita connection nabla_g associated with the Riemannian metric g on the 3-dimensional differentiable manifold M\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "Levi-Civita connection nabla_g associated with the Riemannian metric g on the 3-dimensional differentiable manifold M" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "nabla = g.connection()\n", "print(nabla)\n", "nabla" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The connection coefficients with respect to the spherical orthonormal frame $e$ are" ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "Gam^1_22 = -1/r \n", "Gam^1_33 = -1/r \n", "Gam^2_12 = 1/r \n", "Gam^2_33 = -cos(th)/(r*sin(th)) \n", "Gam^3_13 = 1/r \n", "Gam^3_23 = cos(th)/(r*sin(th)) " ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "nabla.display()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "while those with respect to the natural frame of spherical coordinates (Christoffel symbols) are:" ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "Gam^r_th,th = -r \n", "Gam^r_ph,ph = -r*sin(th)^2 \n", "Gam^th_r,th = 1/r \n", "Gam^th_th,r = 1/r \n", "Gam^th_ph,ph = -cos(th)*sin(th) \n", "Gam^ph_r,ph = 1/r \n", "Gam^ph_th,ph = cos(th)/sin(th) \n", "Gam^ph_ph,r = 1/r \n", "Gam^ph_ph,th = cos(th)/sin(th) " ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "nabla.display(spher.frame())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The covariant derivative of the displacement vector $U$ is" ] }, { "cell_type": "code", "execution_count": 23, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Tensor field nabla_g(U) of type (1,1) on the 3-dimensional differentiable manifold M\n" ] } ], "source": [ "nabU = nabla(U)\n", "print(nabU)" ] }, { "cell_type": "code", "execution_count": 24, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "nabla_g(U) = d(U_1)/dr e_1*e^1 - (U_2(r, th, ph) - d(U_1)/dth)/r e_1*e^2 - (U_3(r, th, ph)*sin(th) - d(U_1)/dph)/(r*sin(th)) e_1*e^3 + d(U_2)/dr e_2*e^1 + (U_1(r, th, ph) + d(U_2)/dth)/r e_2*e^2 - (U_3(r, th, ph)*cos(th) - d(U_2)/dph)/(r*sin(th)) e_2*e^3 + d(U_3)/dr e_3*e^1 + d(U_3)/dth/r e_3*e^2 + (U_2(r, th, ph)*cos(th) + U_1(r, th, ph)*sin(th) + d(U_3)/dph)/(r*sin(th)) e_3*e^3" ] }, "execution_count": 24, "metadata": {}, "output_type": "execute_result" } ], "source": [ "nabU.display()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We convert it to a tensor field of type (0,2) (i.e. a bilinear form) by lowering the upper index with $g$:" ] }, { "cell_type": "code", "execution_count": 25, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Tensor field of type (0,2) on the 3-dimensional differentiable manifold M\n" ] } ], "source": [ "nabU_form = nabU.down(g)\n", "print(nabU_form)" ] }, { "cell_type": "code", "execution_count": 26, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "d(U_1)/dr e^1*e^1 - (U_2(r, th, ph) - d(U_1)/dth)/r e^1*e^2 - (U_3(r, th, ph)*sin(th) - d(U_1)/dph)/(r*sin(th)) e^1*e^3 + d(U_2)/dr e^2*e^1 + (U_1(r, th, ph) + d(U_2)/dth)/r e^2*e^2 - (U_3(r, th, ph)*cos(th) - d(U_2)/dph)/(r*sin(th)) e^2*e^3 + d(U_3)/dr e^3*e^1 + d(U_3)/dth/r e^3*e^2 + (U_2(r, th, ph)*cos(th) + U_1(r, th, ph)*sin(th) + d(U_3)/dph)/(r*sin(th)) e^3*e^3" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ "nabU_form.display()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The **strain tensor** $\\varepsilon$ is defined as the symmetrized part of this tensor:" ] }, { "cell_type": "code", "execution_count": 27, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Field of symmetric bilinear forms on the 3-dimensional differentiable manifold M\n" ] } ], "source": [ "E = nabU_form.symmetrize()\n", "print(E)" ] }, { "cell_type": "code", "execution_count": 28, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "E = d(U_1)/dr e^1*e^1 + 1/2*(r*d(U_2)/dr - U_2(r, th, ph) + d(U_1)/dth)/r e^1*e^2 + 1/2*((r*d(U_3)/dr - U_3(r, th, ph))*sin(th) + d(U_1)/dph)/(r*sin(th)) e^1*e^3 + 1/2*(r*d(U_2)/dr - U_2(r, th, ph) + d(U_1)/dth)/r e^2*e^1 + (U_1(r, th, ph) + d(U_2)/dth)/r e^2*e^2 - 1/2*(U_3(r, th, ph)*cos(th) - sin(th)*d(U_3)/dth - d(U_2)/dph)/(r*sin(th)) e^2*e^3 + 1/2*((r*d(U_3)/dr - U_3(r, th, ph))*sin(th) + d(U_1)/dph)/(r*sin(th)) e^3*e^1 - 1/2*(U_3(r, th, ph)*cos(th) - sin(th)*d(U_3)/dth - d(U_2)/dph)/(r*sin(th)) e^3*e^2 + (U_2(r, th, ph)*cos(th) + U_1(r, th, ph)*sin(th) + d(U_3)/dph)/(r*sin(th)) e^3*e^3" ] }, "execution_count": 28, "metadata": {}, "output_type": "execute_result" } ], "source": [ "E.set_name('E', latex_name=r'\\varepsilon')\n", "E.display()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us display the components of $\\varepsilon$, skipping those that can be deduced by symmetry:" ] }, { "cell_type": "code", "execution_count": 29, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "E_11 = d(U_1)/dr \n", "E_12 = 1/2*(r*d(U_2)/dr - U_2(r, th, ph) + d(U_1)/dth)/r \n", "E_13 = 1/2*((r*d(U_3)/dr - U_3(r, th, ph))*sin(th) + d(U_1)/dph)/(r*sin(th)) \n", "E_22 = (U_1(r, th, ph) + d(U_2)/dth)/r \n", "E_23 = -1/2*(U_3(r, th, ph)*cos(th) - sin(th)*d(U_3)/dth - d(U_2)/dph)/(r*sin(th)) \n", "E_33 = (U_2(r, th, ph)*cos(th) + U_1(r, th, ph)*sin(th) + d(U_3)/dph)/(r*sin(th)) " ] }, "execution_count": 29, "metadata": {}, "output_type": "execute_result" } ], "source": [ "E.display_comp(only_nonredundant=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Stress tensor and Hooke's law\n", "\n", "To form the stress tensor according to Hooke's law, we introduce first the LamÃ© constants:" ] }, { "cell_type": "code", "execution_count": 30, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "ll" ] }, "execution_count": 30, "metadata": {}, "output_type": "execute_result" } ], "source": [ "var('ll', latex_name=r'\\lambda')" ] }, { "cell_type": "code", "execution_count": 31, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "mu" ] }, "execution_count": 31, "metadata": {}, "output_type": "execute_result" } ], "source": [ "var('mu', latex_name=r'\\mu')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The trace (with respect to $g$) of the bilinear form $\\varepsilon$ is obtained by (i) raising the first index (pos=0) by means of $g$ and (ii) by taking the trace of the resulting endomorphism:" ] }, { "cell_type": "code", "execution_count": 32, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Scalar field on the 3-dimensional differentiable manifold M\n" ] } ], "source": [ "trE = E.up(g, pos=0).trace()\n", "print(trE)" ] }, { "cell_type": "code", "execution_count": 33, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "M --> R\n", "(r, th, ph) |--> (U_2(r, th, ph)*cos(th) + (r*d(U_1)/dr + 2*U_1(r, th, ph) + d(U_2)/dth)*sin(th) + d(U_3)/dph)/(r*sin(th))" ] }, "execution_count": 33, "metadata": {}, "output_type": "execute_result" } ], "source": [ "trE.display()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The **stress tensor** $S$ is obtained via Hooke's law for isotropic material:\n", "$$S = \\lambda \\, \\mathrm{tr}\\varepsilon \\; g + 2\\mu \\, \\varepsilon$$" ] }, { "cell_type": "code", "execution_count": 34, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Field of symmetric bilinear forms on the 3-dimensional differentiable manifold M\n" ] } ], "source": [ "S = ll*trE*g + 2*mu*E\n", "print(S)" ] }, { "cell_type": "code", "execution_count": 35, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "S = (ll*U_2(r, th, ph)*cos(th) + ((ll + 2*mu)*r*d(U_1)/dr + 2*ll*U_1(r, th, ph) + ll*d(U_2)/dth)*sin(th) + ll*d(U_3)/dph)/(r*sin(th)) e^1*e^1 + (mu*r*d(U_2)/dr - mu*U_2(r, th, ph) + mu*d(U_1)/dth)/r e^1*e^2 + ((mu*r*d(U_3)/dr - mu*U_3(r, th, ph))*sin(th) + mu*d(U_1)/dph)/(r*sin(th)) e^1*e^3 + (mu*r*d(U_2)/dr - mu*U_2(r, th, ph) + mu*d(U_1)/dth)/r e^2*e^1 + (ll*U_2(r, th, ph)*cos(th) + (ll*r*d(U_1)/dr + 2*(ll + mu)*U_1(r, th, ph) + (ll + 2*mu)*d(U_2)/dth)*sin(th) + ll*d(U_3)/dph)/(r*sin(th)) e^2*e^2 - (mu*U_3(r, th, ph)*cos(th) - mu*sin(th)*d(U_3)/dth - mu*d(U_2)/dph)/(r*sin(th)) e^2*e^3 + ((mu*r*d(U_3)/dr - mu*U_3(r, th, ph))*sin(th) + mu*d(U_1)/dph)/(r*sin(th)) e^3*e^1 - (mu*U_3(r, th, ph)*cos(th) - mu*sin(th)*d(U_3)/dth - mu*d(U_2)/dph)/(r*sin(th)) e^3*e^2 + ((ll + 2*mu)*U_2(r, th, ph)*cos(th) + (ll*r*d(U_1)/dr + 2*(ll + mu)*U_1(r, th, ph) + ll*d(U_2)/dth)*sin(th) + (ll + 2*mu)*d(U_3)/dph)/(r*sin(th)) e^3*e^3" ] }, "execution_count": 35, "metadata": {}, "output_type": "execute_result" } ], "source": [ "S.set_name('S')\n", "S.display()" ] }, { "cell_type": "code", "execution_count": 36, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "S_11 = (ll*U_2(r, th, ph)*cos(th) + ((ll + 2*mu)*r*d(U_1)/dr + 2*ll*U_1(r, th, ph) + ll*d(U_2)/dth)*sin(th) + ll*d(U_3)/dph)/(r*sin(th)) \n", "S_12 = (mu*r*d(U_2)/dr - mu*U_2(r, th, ph) + mu*d(U_1)/dth)/r \n", "S_13 = ((mu*r*d(U_3)/dr - mu*U_3(r, th, ph))*sin(th) + mu*d(U_1)/dph)/(r*sin(th)) \n", "S_22 = (ll*U_2(r, th, ph)*cos(th) + (ll*r*d(U_1)/dr + 2*(ll + mu)*U_1(r, th, ph) + (ll + 2*mu)*d(U_2)/dth)*sin(th) + ll*d(U_3)/dph)/(r*sin(th)) \n", "S_23 = -(mu*U_3(r, th, ph)*cos(th) - mu*sin(th)*d(U_3)/dth - mu*d(U_2)/dph)/(r*sin(th)) \n", "S_33 = ((ll + 2*mu)*U_2(r, th, ph)*cos(th) + (ll*r*d(U_1)/dr + 2*(ll + mu)*U_1(r, th, ph) + ll*d(U_2)/dth)*sin(th) + (ll + 2*mu)*d(U_3)/dph)/(r*sin(th)) " ] }, "execution_count": 36, "metadata": {}, "output_type": "execute_result" } ], "source": [ "S.display_comp(only_nonredundant=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Each component can be accessed individually:" ] }, { "cell_type": "code", "execution_count": 37, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "(mu*r*d(U_2)/dr - mu*U_2(r, th, ph) + mu*d(U_1)/dth)/r" ] }, "execution_count": 37, "metadata": {}, "output_type": "execute_result" } ], "source": [ "S[1,2]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Divergence of the stress tensor\n", "\n", "The divergence of the stress tensor (with respect to $g$) is the 1-form:\n", "$$f_i = \\nabla_j S^j_{\\ \\, i}$$\n", "In a next version of SageManifolds, there will be a function divergence(). For the moment, to evaluate $f$, \n", "we first form the tensor $S^j_{\\ \\, i}$ by raising the first index (pos=0) of $S$ with $g$:" ] }, { "cell_type": "code", "execution_count": 38, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Tensor field of type (1,1) on the 3-dimensional differentiable manifold M\n" ] } ], "source": [ "SU = S.up(g, pos=0)\n", "print(SU)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The divergence is obtained by taking the trace on the first index (0) and the third one (2) of the tensor\n", "$(\\nabla S)^j_{\\ \\, ik} = \\nabla_k S^j_{\\ \\, i}$:" ] }, { "cell_type": "code", "execution_count": 39, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1-form on the 3-dimensional differentiable manifold M\n" ] } ], "source": [ "divS = nabla(SU).trace(0,2)\n", "print(divS)" ] }, { "cell_type": "code", "execution_count": 40, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "f = (((ll + 2*mu)*r^2*d^2(U_1)/dr^2 + 2*(ll + 2*mu)*r*d(U_1)/dr + (ll + mu)*r*d^2(U_2)/drdth - 2*(ll + 2*mu)*U_1(r, th, ph) + mu*d^2(U_1)/dth^2 - (ll + 3*mu)*d(U_2)/dth)*sin(th)^2 + ((ll + mu)*r*d^2(U_3)/drdph + ((ll + mu)*r*d(U_2)/dr - (ll + 3*mu)*U_2(r, th, ph) + mu*d(U_1)/dth)*cos(th) - (ll + 3*mu)*d(U_3)/dph)*sin(th) + mu*d^2(U_1)/dph^2)/(r^2*sin(th)^2) e^1 + ((mu*r^2*d^2(U_2)/dr^2 + (ll + mu)*r*d^2(U_1)/drdth + 2*mu*r*d(U_2)/dr + 2*(ll + 2*mu)*d(U_1)/dth + (ll + 2*mu)*d^2(U_2)/dth^2)*sin(th)^2 - (ll + 3*mu)*cos(th)*d(U_3)/dph - (ll + 2*mu)*U_2(r, th, ph) + ((ll + 2*mu)*cos(th)*d(U_2)/dth + (ll + mu)*d^2(U_3)/dthdph)*sin(th) + mu*d^2(U_2)/dph^2)/(r^2*sin(th)^2) e^2 + ((mu*r^2*d^2(U_3)/dr^2 + 2*mu*r*d(U_3)/dr + mu*d^2(U_3)/dth^2)*sin(th)^2 + (ll + 3*mu)*cos(th)*d(U_2)/dph - mu*U_3(r, th, ph) + ((ll + mu)*r*d^2(U_1)/drdph + mu*cos(th)*d(U_3)/dth + 2*(ll + 2*mu)*d(U_1)/dph + (ll + mu)*d^2(U_2)/dthdph)*sin(th) + (ll + 2*mu)*d^2(U_3)/dph^2)/(r^2*sin(th)^2) e^3" ] }, "execution_count": 40, "metadata": {}, "output_type": "execute_result" } ], "source": [ "divS.set_name('f')\n", "divS.display()" ] }, { "cell_type": "code", "execution_count": 41, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "f_1 = (((ll + 2*mu)*r^2*d^2(U_1)/dr^2 + 2*(ll + 2*mu)*r*d(U_1)/dr + (ll + mu)*r*d^2(U_2)/drdth - 2*(ll + 2*mu)*U_1(r, th, ph) + mu*d^2(U_1)/dth^2 - (ll + 3*mu)*d(U_2)/dth)*sin(th)^2 + ((ll + mu)*r*d^2(U_3)/drdph + ((ll + mu)*r*d(U_2)/dr - (ll + 3*mu)*U_2(r, th, ph) + mu*d(U_1)/dth)*cos(th) - (ll + 3*mu)*d(U_3)/dph)*sin(th) + mu*d^2(U_1)/dph^2)/(r^2*sin(th)^2) \n", "f_2 = ((mu*r^2*d^2(U_2)/dr^2 + (ll + mu)*r*d^2(U_1)/drdth + 2*mu*r*d(U_2)/dr + 2*(ll + 2*mu)*d(U_1)/dth + (ll + 2*mu)*d^2(U_2)/dth^2)*sin(th)^2 - (ll + 3*mu)*cos(th)*d(U_3)/dph - (ll + 2*mu)*U_2(r, th, ph) + ((ll + 2*mu)*cos(th)*d(U_2)/dth + (ll + mu)*d^2(U_3)/dthdph)*sin(th) + mu*d^2(U_2)/dph^2)/(r^2*sin(th)^2) \n", "f_3 = ((mu*r^2*d^2(U_3)/dr^2 + 2*mu*r*d(U_3)/dr + mu*d^2(U_3)/dth^2)*sin(th)^2 + (ll + 3*mu)*cos(th)*d(U_2)/dph - mu*U_3(r, th, ph) + ((ll + mu)*r*d^2(U_1)/drdph + mu*cos(th)*d(U_3)/dth + 2*(ll + 2*mu)*d(U_1)/dph + (ll + mu)*d^2(U_2)/dthdph)*sin(th) + (ll + 2*mu)*d^2(U_3)/dph^2)/(r^2*sin(th)^2) " ] }, "execution_count": 41, "metadata": {}, "output_type": "execute_result" } ], "source": [ "divS.display_comp()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that $f_1$ is quite badly displayed. We get a better view by displaying the components one by one:" ] }, { "cell_type": "code", "execution_count": 42, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "(((ll + 2*mu)*r^2*d^2(U_1)/dr^2 + 2*(ll + 2*mu)*r*d(U_1)/dr + (ll + mu)*r*d^2(U_2)/drdth - 2*(ll + 2*mu)*U_1(r, th, ph) + mu*d^2(U_1)/dth^2 - (ll + 3*mu)*d(U_2)/dth)*sin(th)^2 + ((ll + mu)*r*d^2(U_3)/drdph + ((ll + mu)*r*d(U_2)/dr - (ll + 3*mu)*U_2(r, th, ph) + mu*d(U_1)/dth)*cos(th) - (ll + 3*mu)*d(U_3)/dph)*sin(th) + mu*d^2(U_1)/dph^2)/(r^2*sin(th)^2)" ] }, "execution_count": 42, "metadata": {}, "output_type": "execute_result" } ], "source": [ "divS[1]" ] }, { "cell_type": "code", "execution_count": 43, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "((mu*r^2*d^2(U_2)/dr^2 + (ll + mu)*r*d^2(U_1)/drdth + 2*mu*r*d(U_2)/dr + 2*(ll + 2*mu)*d(U_1)/dth + (ll + 2*mu)*d^2(U_2)/dth^2)*sin(th)^2 - (ll + 3*mu)*cos(th)*d(U_3)/dph - (ll + 2*mu)*U_2(r, th, ph) + ((ll + 2*mu)*cos(th)*d(U_2)/dth + (ll + mu)*d^2(U_3)/dthdph)*sin(th) + mu*d^2(U_2)/dph^2)/(r^2*sin(th)^2)" ] }, "execution_count": 43, "metadata": {}, "output_type": "execute_result" } ], "source": [ "divS[2]" ] }, { "cell_type": "code", "execution_count": 44, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "((mu*r^2*d^2(U_3)/dr^2 + 2*mu*r*d(U_3)/dr + mu*d^2(U_3)/dth^2)*sin(th)^2 + (ll + 3*mu)*cos(th)*d(U_2)/dph - mu*U_3(r, th, ph) + ((ll + mu)*r*d^2(U_1)/drdph + mu*cos(th)*d(U_3)/dth + 2*(ll + 2*mu)*d(U_1)/dph + (ll + mu)*d^2(U_2)/dthdph)*sin(th) + (ll + 2*mu)*d^2(U_3)/dph^2)/(r^2*sin(th)^2)" ] }, "execution_count": 44, "metadata": {}, "output_type": "execute_result" } ], "source": [ "divS[3]" ] } ], "metadata": { "kernelspec": { "display_name": "SageMath 7.5.1", "language": "", "name": "sagemath" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.13" } }, "nbformat": 4, "nbformat_minor": 0 }