*This worksheet illustrates some features of SageManifolds (v0.9) on computations regarding elasticity theory in spherical coordinates.*

First we set up the notebook to display mathematical objects using LaTeX formatting:

In [1]:

```
%display latex
```

We then introduce the Euclidean space as a 3-dimensional differentiable manifold:

In [2]:

```
M = Manifold(3, 'M', start_index=1)
print M
```

We shall make use of spherical coordinates $(r,\theta,\phi)$:

In [3]:

```
spher.<r,th,ph> = M.chart(r'r:(0,+oo) th:(0,pi):\theta ph:(0,2*pi):\phi')
print spher
spher
```

Out[3]:

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.

The natural vector frame of spherical coordinates is

In [4]:

```
spher.frame()
```

Out[4]:

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:

In [5]:

```
to_orthonormal = M.automorphism_field()
to_orthonormal[1,1] = 1
to_orthonormal[2,2] = 1/r
to_orthonormal[3,3] = 1/(r*sin(th))
to_orthonormal.display()
```

Out[5]:

In other words, the change-of-basis matrix is

In [6]:

```
to_orthonormal[:]
```

Out[6]:

We construct the orthonormal frame from the natural frame of spherical coordinates by this change of basis:

In [7]:

```
e = spher.frame().new_frame(to_orthonormal, 'e')
e
```

Out[7]:

In [8]:

```
e[1].display()
```

Out[8]:

In [9]:

```
e[2].display()
```

Out[9]:

In [10]:

```
e[3].display()
```

Out[10]:

At this stage, the default vector frame on $M$ is the first one introduced, namely the natural frame of spherical coordinates:

In [11]:

```
M.default_frame()
```

Out[11]:

Since we prefer the orthonormal frame, we declare

In [12]:

```
M.set_default_frame(e)
```

Then, by default, all vector and tensor fields are displayed with respect to that frame:

In [13]:

```
e[3].display()
```

Out[13]:

To get the same output as in `Out[10]`

, one should specify the frame for display, since this is no longer the default one:

In [14]:

```
e[3].display(spher.frame())
```

Out[14]:

Let define the **displacement vector** $U$ in terms of its components w.r.t. the orthonormal spherical frame:

In [15]:

```
U = M.vector_field(name='U')
U[:] = [function('U_1')(r,th,ph), function('U_2')(r,th,ph), function('U_3')(r,th,ph)]
U.display()
```

Out[15]:

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$):

In [16]:

```
g = M.riemannian_metric('g')
print g
```

Since $e$ is supposed to be an orthonormal frame, we declare that the components of $g$ with respect to it are $\mathrm{diag}(1,1,1)$:

In [17]:

```
g[1,1], g[2,2], g[3,3] = 1, 1, 1
g.display()
```

Out[17]:

The expression of $g$ with respect to the natural frame of spherical coordinates is then

In [18]:

```
g.display(spher.frame())
```

Out[18]:

The covariant derivative operator $\nabla$ is introduced as the (Levi-Civita) connection associated with $g$:

In [19]:

```
nabla = g.connection()
print nabla
nabla
```

Out[19]:

The connection coefficients with respect to the spherical orthonormal frame $e$ are

In [20]:

```
nabla.display()
```

Out[20]:

while those with respect to the natural frame of spherical coordinates (Christoffel symbols) are:

In [21]:

```
nabla.display(spher.frame())
```

Out[21]:

The covariant derivative of the displacement vector $U$ is

In [22]:

```
nabU = nabla(U)
print nabU
```

In [23]:

```
nabU.display()
```

Out[23]:

We convert it to a tensor field of type (0,2) (i.e. a bilinear form) by lowering the upper index with $g$:

In [24]:

```
nabU_form = nabU.down(g)
print nabU_form
```

In [25]:

```
nabU_form.display()
```

Out[25]:

The **strain tensor** $\varepsilon$ is defined as the symmetrized part of this tensor:

In [26]:

```
E = nabU_form.symmetrize()
print E
```

In [27]:

```
E.set_name('E', latex_name=r'\varepsilon')
E.display()
```

Out[27]:

Let us display the components of $\varepsilon$, skipping those that can be deduced by symmetry:

In [28]:

```
E.display_comp(only_nonredundant=True)
```

Out[28]:

To form the stress tensor according to Hooke's law, we introduce first the LamÃ© constants:

In [29]:

```
var('ll', latex_name=r'\lambda')
```

Out[29]:

In [30]:

```
var('mu', latex_name=r'\mu')
```

Out[30]:

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:

In [31]:

```
trE = E.up(g, pos=0).trace()
print trE
```

In [32]:

```
trE.display()
```

Out[32]:

The **stress tensor** $S$ is obtained via Hooke's law for isotropic material:
$$ S = \lambda \, \mathrm{tr}\varepsilon \; g + 2\mu \, \varepsilon$$

In [33]:

```
S = ll*trE*g + 2*mu*E
print S
```

In [34]:

```
S.set_name('S')
S.display()
```

Out[34]:

In [35]:

```
S.display_comp(only_nonredundant=True)
```

Out[35]:

Each component can be accessed individually:

In [36]:

```
S[1,2]
```

Out[36]:

The divergence of the stress tensor (with respect to $g$) is the 1-form:
$$ f_i = \nabla_j S^j_{\ \, i} $$
In a next version of SageManifolds, there will be a function `divergence()`

. For the moment, to evaluate $f$,
we first form the tensor $S^j_{\ \, i}$ by raising the first index (`pos=0`

) of $S$ with $g$:

In [37]:

```
SU = S.up(g, pos=0)
print SU
```

The divergence is obtained by taking the trace on the first index (`0`

) and the third one (`2`

) of the tensor
$(\nabla S)^j_{\ \, ik} = \nabla_k S^j_{\ \, i}$:

In [38]:

```
divS = nabla(SU).trace(0,2)
print divS
```

In [39]:

```
divS.set_name('f')
divS.display()
```

Out[39]:

In [40]:

```
divS.display_comp()
```

Out[40]:

Note that $f_1$ is quite badly displayed. We get a better view by displaying the components one by one:

In [41]:

```
divS[1]
```

Out[41]:

In [42]:

```
divS[2]
```

Out[42]:

In [43]:

```
divS[3]
```

Out[43]:

In [ ]:

```
```