{
"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
}