{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Strain and stress tensors in Cartesian 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_Cartesian.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 Cartesian 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 then introduce the Cartesian coordinates $(x,y,z)$ as a chart $X$ on $M$:"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Chart (M, (x, y, z))\n"
]
},
{
"data": {
"text/html": [
""
],
"text/plain": [
"Chart (M, (x, y, z))"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"X. = M.chart()\n",
"print(X)\n",
"X"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The associated vector frame is"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"Coordinate frame (M, (d/dx,d/dy,d/dz))"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"X.frame()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We shall expand vector and tensor fields not on this frame, which is the default one on $M$:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"Coordinate frame (M, (d/dx,d/dy,d/dz))"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"M.default_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 Cartesian frame:"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"U = U_x(x, y, z) d/dx + U_y(x, y, z) d/dy + U_z(x, y, z) d/dz"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"U = M.vector_field(name='U')\n",
"U[:] = [function('U_x')(x,y,z), function('U_y')(x,y,z), \n",
" function('U_z')(x,y,z)]\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": 8,
"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": [
"We initialize $g$ by declaring that its components with respect to the frame of Cartesian coordinates are\n",
"$\\mathrm{diag}(1,1,1)$:"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"g = dx*dx + dy*dy + dz*dz"
]
},
"execution_count": 9,
"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 covariant derivative operator $\\nabla$ is introduced as the (Levi-Civita) connection associated with $g$: "
]
},
{
"cell_type": "code",
"execution_count": 10,
"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": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"nabla = g.connection()\n",
"print(nabla)\n",
"nabla"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The covariant derivative of the displacement vector $U$ is"
]
},
{
"cell_type": "code",
"execution_count": 11,
"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": 12,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"nabla_g(U) = d(U_x)/dx d/dx*dx + d(U_x)/dy d/dx*dy + d(U_x)/dz d/dx*dz + d(U_y)/dx d/dy*dx + d(U_y)/dy d/dy*dy + d(U_y)/dz d/dy*dz + d(U_z)/dx d/dz*dx + d(U_z)/dy d/dz*dy + d(U_z)/dz d/dz*dz"
]
},
"execution_count": 12,
"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": 13,
"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": 14,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"d(U_x)/dx dx*dx + d(U_x)/dy dx*dy + d(U_x)/dz dx*dz + d(U_y)/dx dy*dx + d(U_y)/dy dy*dy + d(U_y)/dz dy*dz + d(U_z)/dx dz*dx + d(U_z)/dy dz*dy + d(U_z)/dz dz*dz"
]
},
"execution_count": 14,
"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": 15,
"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": 16,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"E = d(U_x)/dx dx*dx + (1/2*d(U_x)/dy + 1/2*d(U_y)/dx) dx*dy + (1/2*d(U_x)/dz + 1/2*d(U_z)/dx) dx*dz + (1/2*d(U_x)/dy + 1/2*d(U_y)/dx) dy*dx + d(U_y)/dy dy*dy + (1/2*d(U_y)/dz + 1/2*d(U_z)/dy) dy*dz + (1/2*d(U_x)/dz + 1/2*d(U_z)/dx) dz*dx + (1/2*d(U_y)/dz + 1/2*d(U_z)/dy) dz*dy + d(U_z)/dz dz*dz"
]
},
"execution_count": 16,
"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": 17,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"E_xx = d(U_x)/dx \n",
"E_xy = 1/2*d(U_x)/dy + 1/2*d(U_y)/dx \n",
"E_xz = 1/2*d(U_x)/dz + 1/2*d(U_z)/dx \n",
"E_yy = d(U_y)/dy \n",
"E_yz = 1/2*d(U_y)/dz + 1/2*d(U_z)/dy \n",
"E_zz = d(U_z)/dz "
]
},
"execution_count": 17,
"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": 18,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"ll"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"var('ll', latex_name=r'\\lambda')"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"mu"
]
},
"execution_count": 19,
"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": 20,
"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": 21,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"M --> R\n",
"(x, y, z) |--> d(U_x)/dx + d(U_y)/dy + d(U_z)/dz"
]
},
"execution_count": 21,
"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": 22,
"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": 23,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"S = ((ll + 2*mu)*d(U_x)/dx + ll*d(U_y)/dy + ll*d(U_z)/dz) dx*dx + (mu*d(U_x)/dy + mu*d(U_y)/dx) dx*dy + (mu*d(U_x)/dz + mu*d(U_z)/dx) dx*dz + (mu*d(U_x)/dy + mu*d(U_y)/dx) dy*dx + (ll*d(U_x)/dx + (ll + 2*mu)*d(U_y)/dy + ll*d(U_z)/dz) dy*dy + (mu*d(U_y)/dz + mu*d(U_z)/dy) dy*dz + (mu*d(U_x)/dz + mu*d(U_z)/dx) dz*dx + (mu*d(U_y)/dz + mu*d(U_z)/dy) dz*dy + (ll*d(U_x)/dx + ll*d(U_y)/dy + (ll + 2*mu)*d(U_z)/dz) dz*dz"
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"S.set_name('S')\n",
"S.display()"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"S_xx = (ll + 2*mu)*d(U_x)/dx + ll*d(U_y)/dy + ll*d(U_z)/dz \n",
"S_xy = mu*d(U_x)/dy + mu*d(U_y)/dx \n",
"S_xz = mu*d(U_x)/dz + mu*d(U_z)/dx \n",
"S_yy = ll*d(U_x)/dx + (ll + 2*mu)*d(U_y)/dy + ll*d(U_z)/dz \n",
"S_yz = mu*d(U_y)/dz + mu*d(U_z)/dy \n",
"S_zz = ll*d(U_x)/dx + ll*d(U_y)/dy + (ll + 2*mu)*d(U_z)/dz "
]
},
"execution_count": 24,
"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": 25,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"mu*d(U_x)/dy + mu*d(U_y)/dx"
]
},
"execution_count": 25,
"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 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": 26,
"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": 27,
"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": 28,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"f = ((ll + 2*mu)*d^2(U_x)/dx^2 + mu*d^2(U_x)/dy^2 + mu*d^2(U_x)/dz^2 + (ll + mu)*d^2(U_y)/dxdy + (ll + mu)*d^2(U_z)/dxdz) dx + ((ll + mu)*d^2(U_x)/dxdy + mu*d^2(U_y)/dx^2 + (ll + 2*mu)*d^2(U_y)/dy^2 + mu*d^2(U_y)/dz^2 + (ll + mu)*d^2(U_z)/dydz) dy + ((ll + mu)*d^2(U_x)/dxdz + (ll + mu)*d^2(U_y)/dydz + mu*d^2(U_z)/dx^2 + mu*d^2(U_z)/dy^2 + (ll + 2*mu)*d^2(U_z)/dz^2) dz"
]
},
"execution_count": 28,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"divS.set_name('f')\n",
"divS.display()"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"f_x = (ll + 2*mu)*d^2(U_x)/dx^2 + mu*d^2(U_x)/dy^2 + mu*d^2(U_x)/dz^2 + (ll + mu)*d^2(U_y)/dxdy + (ll + mu)*d^2(U_z)/dxdz \n",
"f_y = (ll + mu)*d^2(U_x)/dxdy + mu*d^2(U_y)/dx^2 + (ll + 2*mu)*d^2(U_y)/dy^2 + mu*d^2(U_y)/dz^2 + (ll + mu)*d^2(U_z)/dydz \n",
"f_z = (ll + mu)*d^2(U_x)/dxdz + (ll + mu)*d^2(U_y)/dydz + mu*d^2(U_z)/dx^2 + mu*d^2(U_z)/dy^2 + (ll + 2*mu)*d^2(U_z)/dz^2 "
]
},
"execution_count": 29,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"divS.display_comp()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Displaying the components one by one:"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"(ll + 2*mu)*d^2(U_x)/dx^2 + mu*d^2(U_x)/dy^2 + mu*d^2(U_x)/dz^2 + (ll + mu)*d^2(U_y)/dxdy + (ll + mu)*d^2(U_z)/dxdz"
]
},
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"divS[1]"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"(ll + mu)*d^2(U_x)/dxdy + mu*d^2(U_y)/dx^2 + (ll + 2*mu)*d^2(U_y)/dy^2 + mu*d^2(U_y)/dz^2 + (ll + mu)*d^2(U_z)/dydz"
]
},
"execution_count": 31,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"divS[2]"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"(ll + mu)*d^2(U_x)/dxdz + (ll + mu)*d^2(U_y)/dydz + mu*d^2(U_z)/dx^2 + mu*d^2(U_z)/dy^2 + (ll + 2*mu)*d^2(U_z)/dz^2"
]
},
"execution_count": 32,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"divS[3]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
}
],
"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
}