Welcome to Continuum Mechanics’s documentation!

Continuum Mechanics

Utilities for doing calculations in continuum mechanics.

• Free software: MIT license.
• Documentation: https://continuum-mechanics.readthedocs.io.

Features

• Make long vector calculus calculations simple.
• Support different coordinate systems.
• Visualize entities such as second rank tensors.

Citation

To cite continuum_mechanics in publications use

Nicolás Guarín-Zapata. (2020). nicoguaro/continuum_mechanics: Version 0.2.1 (Version v0.2.1). Zenodo. http://doi.org/10.5281/zenodo.4029448

A BibTeX entry for LaTeX users is

@software{continuum_mechanics,
 title = {continuum_mechanics: Continuum Mechanics calculations in Python},
 version = {0.2.1},
 author = {Guarín-Zapata, Nicolás},
 year = 2020,
 keywords = {Python, Finite elements, Scientific computing, Computational mechanics},
 abstract = {`continuum_mechanics` is a Python package built on top of SymPy to aid
    with calculations in Continuum Mechanics that are commonly lengthy and
    tedious if done by hand. It also provides visualization capabilities for
    second-order tensors such as Mohr's circle to help in stress analyses.},
 url = {https://github.com/nicoguaro/continuum_mechanics},
 doi = {http://doi.org/10.5281/zenodo.4029448}
}

Credits

This package was created with Cookiecutter and the audreyr/cookiecutter-pypackage project template.

Installation

Stable release

To install Continuum Mechanics, run this command in your terminal:

$ pip install continuum_mechanics

This is the preferred method to install Continuum Mechanics, as it will always install the most recent stable release.

If you don’t have pip installed, this Python installation guide can guide you through the process.

From sources

The sources for Continuum Mechanics can be downloaded from the Github repo.

You can either clone the public repository:

$ git clone git://github.com/nicoguaro/continuum_mechanics

Or download the tarball:

$ curl  -OL https://github.com/nicoguaro/continuum_mechanics/tarball/master

Once you have a copy of the source, you can install it with:

$ python setup.py install

Usage

continuum_mechanics relies on SymPy for symbolic calculations. To use Continuum Mechanics in a project you can import it as:

import continuum_mechanics

Vector operators in Cartesian coordinates

continuum_mechanics support major vector operators such as:

• gradient of a scalar function;
• divergence of a vector function;
• curl of a vector function;
• gradient of a vector function;
• divergence of a tensor;
• Laplace operator of a scalar function;
• Laplace operator of a vector function; and
• Biharmonic operator of a scalar function.

All these operators are in the module vector.

from sympy import *
from continuum_mechanics import vector

By default, Cartesian coordinates are given by \(x\), \(y\) and \(z\). If these coordinates are used there is not necessary to specify them when calling the vector operators

init_printing()
x, y, z = symbols("x y z")

Following, we have some examples of vector operators in Cartesian coordinates.

Gradient of a scalar function

The gradient takes as input a scalar and returns a vector, represented by a 3 by 1 matrix.

f = 2*x + 3*y**2 - sin(z)
f

\[2 x + 3 y^{2} - \sin{\left (z \right )}\]

vector.grad(f)

\[\begin{split}\left[\begin{matrix}2\\6 y\\- \cos{\left (z \right )}\end{matrix}\right]\end{split}\]

Divergence of a vector function

The divergence takes as input a vector (represented by a 3 by 1 matrix) and returns a scalar.

vector.div(Matrix([x, y, z]))

\[3\]

vector.div(Matrix([
    x**2 + y*z,
    y**2 + x*z,
    z**2 + x*y]))

\[2 x + 2 y + 2 z\]

Divergence of a tensor function

The divergence of a tensor (represented by a 3 by 3 matrix) returns a vector.

Axx, Axy, Axz = symbols("A_xx A_xy A_xz", cls=Function)
Ayx, Ayy, Ayz = symbols("A_yx A_yy A_yz", cls=Function)
Azx, Azy, Azz = symbols("A_zx A_zy A_zz", cls=Function)

tensor = Matrix([
    [Axx(x, y, z), Axy(x, y, z), Axz(x, y, z)],
    [Ayx(x, y, z), Ayy(x, y, z), Ayz(x, y, z)],
    [Azx(x, y, z), Azy(x, y, z), Azz(x, y, z)]])
tensor

\[\begin{split}\begin{bmatrix} A_{xx}(x,y,z) & A_{xy}(x,y,z) & A_{xz}(x,y,z)\\ A_{yx}(x,y,z) & A_{yy}(x,y,z) & A_{yz}(x,y,z)\\ A_{zx}(x,y,z) & A_{zy}(x,y,z) & A_{zz}(x,y,z) \end{bmatrix}\end{split}\]

vector.div_tensor(tensor)

\[\begin{split}\left[\begin{matrix} \frac{\partial}{\partial x} A_{xx}(x, y, z) + \frac{\partial}{\partial y} A_{xy}(x, y, z) + \frac{\partial}{\partial z} A_{xz}(x, y, z)\\ \frac{\partial}{\partial x} A_{yx}(x, y, z) + \frac{\partial}{\partial y} A_{yy}(x, y, z) + \frac{\partial}{\partial z} A_{yz}(x, y, z)\\ \frac{\partial}{\partial x} A_{zx}(x, y, z) + \frac{\partial}{\partial y} A_{zy}(x, y, z) + \frac{\partial}{\partial z} A_{zz}(x, y, z) \end{matrix}\right]\end{split}\]

Curl of a vector function

Let us check the identity

\[\nabla \times \nabla f(x, y, z) = 0\, .\]

fun = symbols("fun", cls=Function)
vector.curl(vector.grad(fun(x, y, z)))

\[\begin{split}\left[\begin{matrix}0\\0\\0\end{matrix}\right]\end{split}\]

Visualization of tensors

from sympy import Matrix
from continuum_mechanics.visualization import mohr2d, mohr3d, traction_circle

Visualization in 2D

First, let us visualize the tensor

\[\begin{split}\begin{bmatrix} 1 &0\\ 0 &-1 \end{bmatrix}\, .\end{split}\]

mohr2d(Matrix([
  [1,0],
  [0,-1]]))

From the Mohr circle, we can see that the principal directions are given at \(0\) and \(\pi/2\) radians. This can be more easily visualized using the traction circle, where normal vectors are presented in light gray and the traction vectors are presented in colors.

traction_circle(Matrix([
  [1,0],
  [0,-1]]))

Now, let us visualize

\[\begin{split}\begin{bmatrix} 1 &3\\ 3 &-5 \end{bmatrix}\, .\end{split}\]

mohr2d(Matrix([
  [1, 3],
  [3, -5]]))

traction_circle(Matrix([
  [1, 3],
  [3, -5]]))

Now, let us try it with an asymmetric tensor

\[\begin{split}\begin{bmatrix} 1 &2\\ 0 &3 \end{bmatrix}\, .\end{split}\]

mohr2d(Matrix([
  [1, 2],
  [0, 3]]))

traction_circle(Matrix([
  [1, 2],
  [0, 3]]))

Mohr Circle in 3D

Let us visualize the tensor

\[\begin{split}\begin{bmatrix} 1 &2 &4\\ 2 &2 &1\\ 4 &1 &3 \end{bmatrix}\, .\end{split}\]

mohr3d(Matrix([
    [1, 2, 4],
    [2, 2, 1],
    [4, 1, 3]]))

Now, let us visualize the tensor

\[\begin{split}\begin{bmatrix} 1 &0 &0\\ 0 &2 &0\\ 0 &0 &3 \end{bmatrix}\, .\end{split}\]

mohr3d(Matrix([
    [1, 0, 0],
    [0, 2, 0],
    [0, 0, 3]]))

Elasticity tensor visualization

Let us consider β-brass that is a cubic material and has the following material properties in Voigt notation:

\[C_{11} = 52\text{ GPa},\quad C_{12} = 27.5\text{ GPa},\quad C_{44} = 173\text{ GPa}.\]

C11 = 52e9
C12 = 27.5e9
C44 = 173e9
rho = 7600
C = np.zeros((6, 6))
C[0:3, 0:3] = np.array([[C11, C12, C12],
                        [C12, C11, C12],
                        [C12, C12, C11]])
C[3:6, 3:6] = np.diag([C44, C44, C44])

One way to visualize a stiffness tensor is to use the Christofel equation

\[\det(\Gamma_{ij} - v_p^2\delta_{ij}) = 0\, ,\]

where \(\Gamma_{ij}\) is the Christofel stiffness and depends on the material properties (\(c_{ijkl}\)) and unit vectors (:math`n_i`):

This provides the eigenvalues that represent the phase speed for three propagation modes in the material.

V1, V2, V3, phi_vec, theta_vec = christofel_eig(C, 100, 100)
V1 = np.sqrt(V1/rho)
V2 = np.sqrt(V2/rho)
V3 = np.sqrt(V3/rho)

Phase speed for the first quasi-transverse mode.

Phase speed for the second quasi-transverse mode.

Phase speed for the quasi-longitudinal mode.

Continuum Mechanics’ tutorials

Here we show how to use continuum_mechanics’ capabilities through tutorials.

Introduction to orthogonal coordinates

In \(\mathbb{R}^3\), we can think that each point is given by the intersection of three surfaces. Thus, we have three families of curved surfaces that intersect each other at right angles. These surfaces are orthogonal locally, but not (necessarily) globally, and are defined by

\[u_1 = f_1(x, y, z)\, ,\quad u_2 = f_2(x, y, z)\, ,\quad u_3=f_3(x, y, z) \, .\]

These functions should be invertible, at least locally, and we can also write

\[x = x(u_1, u_2, u_3)\, ,\quad y = y(u_1, u_2, u_3)\, ,\quad z = z(u_1, u_2, u_3)\, ,\]

where \(x, y, z\) are the usual Cartesian coordinates. The curve defined by the intersection of two of the surfaces gives us one of the coordinate curves.

Scale factors

Since we are interested in how these surface intersect each other locally, we want to express differential vectors in terms of the coordinates. Thus, the differential for the position vector (\(\mathbf{r}\)) is given by

\[\mathrm{d}\mathbf{r} = \frac{\partial\mathbf{r}}{\partial u_1}\mathrm{d}u_1 + \frac{\partial\mathbf{r}}{\partial u_2}\mathrm{d}u_2 + \frac{\partial\mathbf{r}}{\partial u_3}\mathrm{d}u_3\, ,\]

or

\[\mathrm{d}\mathbf{r} = \sum_{i=1}^3 \frac{\partial\mathbf{r}}{\partial u_i}\mathrm{d}u_i\, .\]

The factor \(\partial \mathbf{r}/\partial u_i\) is a non-unitary vector that takes into account the variation of \(\mathbf{r}\) in the direction of \(u_i\), and is then tangent to the coordinate curve \(u_i\). We can define a normalized basis \(\hat{\mathbf{e}}_i\) using

\[\frac{\partial\mathbf{r}}{\partial u_i} = h_i \hat{\mathbf{e}}_i\, .\]

The coefficients \(h_i\) are functions of \(u_i\) and we call them scale factors. They are really important since they allow us to measure distances while we move along our coordinates. We would need them to define vector operators in orthogonal coordinates. When the coordinates are not orthogonal we would need to use the metric tensor, but we are going to restrict ourselves to orthogonal systems.

Hence, we have the following

\[\begin{split}\begin{align} &h_i = \left|\frac{\partial\mathbf{r}}{\partial u_i}\right|\, ,\\ &\hat{\mathbf{e}}_i = \frac{1}{h_i} \frac{\partial \mathbf{r}}{\partial u_i}\, . \end{align}\end{split}\]

Curvilinear coordinates available

The following coordinate systems are available:

• Cartesian;
• Cylindrical;
• Spherical;
• Parabolic cylindrical;
• Parabolic;
• Paraboloidal;
• Elliptic cylindrical;
• Oblate spheroidal;
• Prolate spheroidal;
• Ellipsoidal;
• Bipolar cylindrical;
• Toroidal;
• Bispherical; and
• Conical.

To obtain the transformation for a given coordinate system we can use the function vector.transform_coords().

import sympy as sym
from continuum_mechanics import vector

First, we define the variables for the coordinates \((u, v, w)\).

sym.init_printing()
u, v, w = sym.symbols("u v w")

And, we compute the coordinates for the parabolic system using vector.transform_coords(). The first parameter is a string defining the coordinate system and the second is a tuple with the coordinates.

vector.transform_coords("parabolic", (u, v, w))

\[\displaystyle \left( u v \cos{\left(w \right)}, \ u v \sin{\left(w \right)}, \ \frac{u^{2}}{2} - \frac{v^{2}}{2}\right)\]

The scale factors for the coordinate systems mentioned above are availabe. We can compute them for bipolar cylindrical coordinates. The coordinates are defined by

\[\begin{split}\begin{align} &x = a \frac{\sinh\tau}{\cosh\tau - \cos\sigma}\, ,\\ &y = a \frac{\sin\sigma}{\cosh\tau - \cos\sigma}\, ,\\ &z = z\, , \end{align}\end{split}\]

and have the following scale factors

\[h_\sigma = h_\tau = \frac{a}{\cosh\tau - \cos\sigma}\, ,\]

and \(h_z = 1\).

sigma, tau, z, a = sym.symbols("sigma tau z a")
z = sym.symbols("z")
scale = vector.scale_coeff_coords("bipolar_cylindrical", (sigma, tau, z), a=a)
scale

\[\displaystyle \left( \frac{a}{- \cos{\left(\sigma \right)} + \cosh{\left(\tau \right)}}, \ \frac{a}{- \cos{\left(\sigma \right)} + \cosh{\left(\tau \right)}}, \ 1\right)\]

Finally, we can compute vector operators for different coordinates.

The Laplace operator for the bipolar cylindrical system is given by

\[\nabla^2 \phi = \frac{1}{a^2} \left( \cosh \tau - \cos\sigma \right)^{2} \left( \frac{\partial^2 \phi}{\partial \sigma^2} + \frac{\partial^2 \phi}{\partial \tau^2} \right) + \frac{\partial^2 \phi}{\partial z^2}\, ,\]

and we can compute it using the function vector.lap(). For this function, the first parameter is the expression that we want to compute the Laplacian for, the second parameter is a tuple with the coordinates and the third parameter is a tuple with the scale factors.

phi = sym.symbols("phi", cls=sym.Function)
lap = vector.lap(phi(sigma, tau, z), coords=(sigma, tau, z), h_vec=scale)
sym.simplify(lap)

\[\displaystyle \frac{a^{2} \frac{\partial^{2}}{\partial z^{2}} \phi{\left(\sigma,\tau,z \right)} + \left(\cos{\left(\sigma \right)} - \cosh{\left(\tau \right)}\right)^{2} \frac{\partial^{2}}{\partial \sigma^{2}} \phi{\left(\sigma,\tau,z \right)} + \left(\cos{\left(\sigma \right)} - \cosh{\left(\tau \right)}\right)^{2} \frac{\partial^{2}}{\partial \tau^{2}} \phi{\left(\sigma,\tau,z \right)}}{a^{2}}\]

Point source on top of a halfspace

import numpy as np
from sympy import init_printing, symbols, lambdify, S, simplify
from sympy import pi, Matrix, sqrt, oo
from continuum_mechanics.solids import sym_grad, strain_stress

import matplotlib.pyplot as plt
from matplotlib import colors

The following snippet allows to format the graphs.

repo = "https://raw.githubusercontent.com/nicoguaro/matplotlib_styles/master"
style = repo + "/styles/minimalist.mplstyle"
plt.style.use(style)

x, y, z, r, E, nu, Fx, Fy, Fz = symbols('x y z r E nu F_x F_y F_z')

The components of the displacement vector are given by [LANDAU]_

\[\begin{split}\begin{align} &u_x = \frac{(1 + \nu)}{2 \pi E} \left\{\left[\frac{xz}{r^3} - \frac{(1 - 2\nu)x}{r(r + z)}\right]F_z + \frac{2(1 - \nu)r +z}{r(r + z)}F_x +\frac{[2r(\nu r + z) + z^2]x}{r^3(r + z)^2}(xF_x + y F_y)\right\}\, ,\\ &u_y = \frac{(1 + \nu)}{2 \pi E} \left\{\left[\frac{yz}{r^3} - \frac{(1 - 2\nu)y}{r(r + z)}\right]F_z + \frac{2(1 - \nu)r +z}{r(r + z)}F_y +\frac{[2r(\nu r + z) + z^2]y}{r^3(r + z)^2}(xF_x + y F_y)\right\}\, ,\\ &u_z = \frac{(1 + \nu)}{2 \pi E} \left\{\left[\frac{2(1 - \nu)}{r} - \frac{z^2}{r^3}\right]F_z +\left[\frac{1 - 2\nu}{r(r + z)} + \frac{z}{r^3}\right](xF_x + y F_y)\right\}\, , \end{align}\end{split}\]

with \(r = \sqrt{x^2 + y^2 + z^2}\).

ux = (1+nu)/(2*pi*E)*((x*z/r**3 - (1-2*nu)*x/(r*(r + z)))*Fz +
                       (2*(1 - nu)*r + z)/(r*(r + z))*Fx +
                       ((2*r*(nu*r + z) + z**2)*x)/(r**3*(r + z)**2)*(x*Fx + y*Fy))
ux

\[\frac{1}{2 \pi E} \left(\nu + 1\right) \left(\frac{F_{x}}{r \left(r + z\right)} \left(r \left(- 2 \nu + 2\right) + z\right) + F_{z} \left(- \frac{x \left(- 2 \nu + 1\right)}{r \left(r + z\right)} + \frac{x z}{r^{3}}\right) + \frac{x}{r^{3} \left(r + z\right)^{2}} \left(F_{x} x + F_{y} y\right) \left(2 r \left(\nu r + z\right) + z^{2}\right)\right)\]

uy = (1+nu)/(2*pi*E)*((y*z/r**3 - (1-2*nu)*y/(r*(r + z)))*Fz +
                       (2*(1 - nu)*r + z)/(r*(r + z))*Fy +
                       ((2*r*(nu*r + z) + z**2)*y)/(r**3*(r + z)**2)*(x*Fx + y*Fy))
uy

\[\frac{1}{2 \pi E} \left(\nu + 1\right) \left(\frac{F_{y}}{r \left(r + z\right)} \left(r \left(- 2 \nu + 2\right) + z\right) + F_{z} \left(- \frac{y \left(- 2 \nu + 1\right)}{r \left(r + z\right)} + \frac{y z}{r^{3}}\right) + \frac{y}{r^{3} \left(r + z\right)^{2}} \left(F_{x} x + F_{y} y\right) \left(2 r \left(\nu r + z\right) + z^{2}\right)\right)\]

uz = (1+nu)/(2*pi*E)*((2*(1 - nu)/r + z**2/r**3)*Fz +
                       ((1 - 2*nu)/(r*(r + z)) + z/r**3)*(x*Fx + y*Fy))
uz

\[\frac{1}{2 \pi E} \left(\nu + 1\right) \left(F_{z} \left(\frac{1}{r} \left(- 2 \nu + 2\right) + \frac{z^{2}}{r^{3}}\right) + \left(F_{x} x + F_{y} y\right) \left(\frac{- 2 \nu + 1}{r \left(r + z\right)} + \frac{z}{r^{3}}\right)\right)\]

Withouth loss of generality we can assume that \(F_y=0\), this is equivalent a rotate the axes until the force is in the plane \(y=0\).

ux = ux.subs(Fy, 0)
ux

\[\frac{1}{2 \pi E} \left(\nu + 1\right) \left(\frac{Fx}{r \left(r + z\right)} \left(r \left(- 2 \nu + 2\right) + z\right) + \frac{Fx x^{2}}{r^{3} \left(r + z\right)^{2}} \left(2 r \left(\nu r + z\right) + z^{2}\right) + Fz \left(- \frac{x \left(- 2 \nu + 1\right)}{r \left(r + z\right)} + \frac{x z}{r^{3}}\right)\right)\]

uy = ux.subs(Fy, 0)
uy

\[\frac{1}{2 \pi E} \left(\nu + 1\right) \left(\frac{Fx}{r \left(r + z\right)} \left(r \left(- 2 \nu + 2\right) + z\right) + \frac{Fx x^{2}}{r^{3} \left(r + z\right)^{2}} \left(2 r \left(\nu r + z\right) + z^{2}\right) + Fz \left(- \frac{x \left(- 2 \nu + 1\right)}{r \left(r + z\right)} + \frac{x z}{r^{3}}\right)\right)\]

uz = uz.subs(Fy, 0)
uz

\[\frac{1}{2 \pi E} \left(\nu + 1\right) \left(Fx x \left(\frac{- 2 \nu + 1}{r \left(r + z\right)} + \frac{z}{r^{3}}\right) + Fz \left(\frac{1}{r} \left(- 2 \nu + 2\right) + \frac{z^{2}}{r^{3}}\right)\right)\]

The displacement vector is then

u = Matrix([ux, uy, uz]).subs(r, sqrt(x**2 + y**2 + z**2))

We can check that the displacement vanish when \(x,y,z \rightarrow \infty\)

u.limit(x, oo)

\[\begin{split}\left[\begin{matrix}0\\0\\0\end{matrix}\right]\end{split}\]

u.limit(y, oo)

\[\begin{split}\left[\begin{matrix}0\\0\\0\end{matrix}\right]\end{split}\]

u.limit(z, oo)

\[\begin{split}\left[\begin{matrix}0\\0\\0\end{matrix}\right]\end{split}\]

We can compute the strain and stress tensors using the symmetric gradient (vector.sym_grad()) and strain-to-stress (solids.strain_stress()) functions.

lamda = E*nu/((1 + nu)*(1 - 2*nu))
mu = E/(2*(1 - nu))
strain = sym_grad(u)
stress = strain_stress(strain, [lamda, mu])

The expressions for strains and stresses are lengthy and difficult to work with. Nevertheless, we can work with them. For example, we can evaluate the stress tensor at a point \(\mathbf{x} = (1, 0, 1)\) for a vertical load and a Poisson coefficient \(\nu = 1/4\).

simplify(stress.subs({x: 1, y: 0, z:1, nu: S(1)/4, Fx: 0}))

\[\begin{split}\left[\begin{matrix}- \frac{F_{z} \left(73 \sqrt{2} + 108\right)}{48 \pi \left(7 + 5 \sqrt{2}\right)} & - \frac{5 F_{z} \left(2 \sqrt{2} + 3\right)}{24 \pi \left(7 + 5 \sqrt{2}\right)} & - \frac{5 F_{z} \left(4 + 3 \sqrt{2}\right)}{16 \pi \left(2 \sqrt{2} + 3\right)}\\- \frac{5 F_{z} \left(2 \sqrt{2} + 3\right)}{24 \pi \left(7 + 5 \sqrt{2}\right)} & - \frac{F_{z} \left(11 \sqrt{2} + 16\right)}{16 \pi \left(7 + 5 \sqrt{2}\right)} & 0\\- \frac{5 F_{z} \left(4 + 3 \sqrt{2}\right)}{16 \pi \left(2 \sqrt{2} + 3\right)} & 0 & - \frac{F_{z} \left(103 \sqrt{2} + 148\right)}{48 \pi \left(7 + 5 \sqrt{2}\right)}\end{matrix}\right]\end{split}\]

Visualization of the fields

Since it is difficult to handle these lengthy expressions we can visualize them. For that, we define a grid where to evaluate the expressions,

\[(x, z) \in [-2, 2]\times[0, 5]\, ,\]

in this case.

x_vec, z_vec = np.mgrid[-2:2:100j, 0:5:100j]

We can use lambdify() to turn the SymPy expressions to evaluatable functions.

def field_plot(expr, x_vec, y_vec, z_vec, E_val, nu_val, Fx_val, Fz_val, title=''):
    """Plot the field"""

    # Lambdify the function
    expr_fun = lambdify((x, y, z, E, nu, Fx, Fz), expr, "numpy")
    expr_vec = expr_fun(x_vec, y_vec, z_vec, E_val, nu_val, Fx_val, Fz_val)

    # Determine extrema
    vmin = np.min(expr_vec)
    vmax = np.max(expr_vec)
    print("Minimum value in the domain: {:g}".format(vmin))
    print("Maximum value in the domain: {:g}".format(vmax))
    vmax = max(np.abs(vmax), np.abs(vmin))

    # Plotting
    fig = plt.gcf()
    levels = np.logspace(-1, np.log10(vmax), 10)
    levels = np.hstack((-levels[-1::-1], [0], levels))
    cbar_ticks = ["{:.2g}".format(level) for level in levels]
    cont = plt.contourf(x_vec, z_vec, expr_vec, levels=levels,
                        cmap="RdYlBu_r", norm=colors.SymLogNorm(0.1))
    cbar = fig.colorbar(cont, ticks=levels[::2])
    cbar.ax.set_yticklabels(cbar_ticks[::2])
    plt.axis("image")
    plt.gca().invert_yaxis()
    plt.xlabel(r"$x$")
    plt.ylabel(r"$z$")
    plt.title(title)
    return cont

Displacements

plt.figure()
field_plot(u.norm(), x_vec, 0, z_vec, 1.0, 0.3, 0.0, 1.0)
plt.show()

Minimum value in the domain: 0.0881197
Maximum value in the domain: 15.4645

plt.figure()
field_plot(u[0], x_vec, 0, z_vec, 1.0, 0.3, 0.0, 1.0)
plt.show()

Minimum value in the domain: -4.09665
Maximum value in the domain: 4.09665

plt.figure()
field_plot(u[2], x_vec, 0, z_vec, 1.0, 0.3, 0.0, 1.0)
plt.show()

Minimum value in the domain: 0.0869101
Maximum value in the domain: 14.3383

Stresses

We can plot the components of stress

for row in range(0, 3):
    for col in range(row, 3):
        plt.figure()
        field_plot(stress[row,col], x_vec, 0, z_vec, 1.0, 0.3, 0.0, 1.0,
                   title=r"$\sigma_{%i%i}$"%(row+1, col+1))

plt.show()

Minimum value in the domain: -41.4274
Maximum value in the domain: 406.682

Minimum value in the domain: -12.0021
Maximum value in the domain: 144.846

Minimum value in the domain: -95.9472
Maximum value in the domain: 95.9472

Minimum value in the domain: -59.0538
Maximum value in the domain: 116.991

Minimum value in the domain: -506.96
Maximum value in the domain: 506.96

Minimum value in the domain: -243.272
Maximum value in the domain: 116.991

Stress invariants

We can also plot the invariants of the stress tensor

I1 = S(1)/3 * stress.trace()
I2 = S(1)/2 * (stress.trace()**2 + (stress**2).trace())
I3 = stress.det()
Mises = sqrt(((stress[0,0] - stress[1,1])**2 + (stress[1,1] - stress[2,2])**2 +
             (stress[2,2] - stress[0,0])**2 +
             6*(stress[0,1]**2 + stress[1,2]**2 + stress[0,2]**2))/2)

plt.figure()
field_plot(I1, x_vec, 0, z_vec, 1.0, 0.3, 0.0, 1.0)
plt.show()

Minimum value in the domain: -107.797
Maximum value in the domain: 213.555

plt.figure()
field_plot(I2, x_vec, 0, z_vec, 1.0, 0.3, 0.0, 1.0)
plt.show()

Minimum value in the domain: 0.000977492
Maximum value in the domain: 579596

plt.figure()
field_plot(I3, x_vec, 0, z_vec, 1.0, 0.3, 0.0, 1.0)
plt.show()

Minimum value in the domain: -1.01409e+08
Maximum value in the domain: 419218

plt.figure()
field_plot(Mises, x_vec, 0, z_vec, 1.0, 0.3, 0.0, 1.0)
plt.show()

Minimum value in the domain: 0.0274784
Maximum value in the domain: 958.065

References

[LANDAU]

Landau, L. D., Kosevich, A. M., Pitaevskii, L. P., & Lifshitz, E. M. (1986). Theory of elasticity.

Dispersion relations in a micropolar medium

We are interested in computing the dispersion relations in a homogeneous micropolar solid.

Wave propagation in micropolar solids

The equations of motion for a micropolar solid are given by [NOWACKI], [GUARIN]

\[\begin{split}\begin{align} &c_1^2 \nabla\nabla\cdot\mathbf{u}- c_2^2\nabla\times\nabla\times\mathbf{u} + K^2\nabla\times\boldsymbol{\theta} = -\omega^2 \mathbf{u} \, ,\\ &c_3^2 \nabla\nabla\cdot\boldsymbol{\theta} - c_4^2\nabla\times\nabla\times\boldsymbol{\theta} + Q^2\nabla\times\mathbf{u} - 2Q^2\boldsymbol{\theta} = -\omega^2 \boldsymbol{\theta} \, \end{align}\end{split}\]

where \(\mathbf{u}\) is the displacement vector and \(\boldsymbol{\theta}\) is the microrrotations vector, and where: \(c_1\) represents the phase/group speed for the longitudinal wave (\(P\)) that is non-dispersive as in the classical case, \(c_2\) represents the high-frequency limit phase/group speed for a transverse wave (\(S\)) that is dispersive unlike the classical counterpart, \(c_3\) represents the high-frequency limit phase/group speed for a longitudinal-rotational wave (\(LR\)) with a corkscrew-like motion that is dispersive and does not have a classical counterpart, \(c_4\) represents the high-frequency limit phase/group speed for a transverse-rotational wave (\(TR\)) that is dispersive and does not have a classical counterpart, \(Q\) represents the cut-off frequency for rotational waves appearance, and \(K\) quantifies the difference between the low-frequency and high-frequency phase/group speed for the S-wave. These parameters are defined by:

\[\begin{split}\begin{align} c_1^2 = \frac{\lambda +2\mu}{\rho},\quad &c_3^2 =\frac{\beta + 2\eta}{J},\\ c_2^2 = \frac{\mu +\alpha}{\rho},\quad &c_4^2 =\frac{\eta + \varepsilon}{J},\\ Q^2= \frac{2\alpha}{J},\quad &K^2 =\frac{2\alpha}{\rho} \, , \end{align}\end{split}\]

Dispersion relations

To identify types of propagating waves that can arise in the micropolar medium it is convenient to expand the displacement and rotation vectors in terms of scalar and vector potentials

\[\begin{split}\begin{align} \mathbf{u} &= \nabla \phi + \nabla\times\boldsymbol{\Gamma}\, ,\\ \boldsymbol{\theta} &= \nabla \tau + \nabla\times\mathbf{E}\, , \end{align}\end{split}\]

subject to the conditions:

\[\begin{split}\begin{align} &\nabla\cdot\boldsymbol{\Gamma} = 0\\ &\nabla\cdot\mathbf{E} = 0\, . \end{align}\end{split}\]

Using the above in the displacements equations of motion yields the following equations, after some manipulations

\[\begin{split}\begin{align} c_1^2 \nabla^2 \phi &= \frac{\partial^2 \phi}{\partial t^2}\, ,\\ c_3^2 \nabla^2 \tau - 2Q^2\tau &= \frac{\partial^2 \tau}{\partial t^2}\, ,\\ \begin{bmatrix} c_2^2 \nabla^2 &K^2\nabla\times\, ,\\ Q^2\nabla\times &c_4^2\nabla^2 - 2Q^2 \end{bmatrix} \begin{Bmatrix} \boldsymbol{\Gamma}\\ \mathbf{E}\end{Bmatrix} &= \frac{\partial^2}{\partial t^2} \begin{Bmatrix} \boldsymbol{\Gamma}\\ \mathbf{E}\end{Bmatrix} \, , \end{align}\end{split}\]

where we can see that the equations for the scalar potentials are uncoupled, while the ones for the vector potentials are coupled.

Writing the vector potentials as plane waves of amplitude \(\mathbf{A}\) and \(\mathbf{B}\), wave number \(\kappa\) and circular frequency \(\omega\) that propagate along the (x) axis,

\[\begin{split}\begin{align} \boldsymbol{\Gamma} &= \mathbf{A}\exp(i\kappa x - i\omega t)\\ \mathbf{E} &= \mathbf{B}\exp(i\kappa x - i\omega t)\, . \end{align}\end{split}\]

We can do these calculations using some the functions available functions in the package.

from sympy import Matrix, diff, symbols, exp, I, sqrt
from sympy import simplify, expand, solve, limit
from sympy import init_printing, pprint, factor
from continuum_mechanics.vector import lap_vec, curl, div

A1, A2, A3, B1, B2, B3 = symbols("A1 A2 A3 B1 B2 B3")
kappa, omega, t, x = symbols("kappa omega t x")
c1, c2, c3, c4, K, Q = symbols("c1 c2 c3 c4 K Q", positive=True)

We define the vector potentials \(\boldsymbol{\Gamma}\) and \(\mathbf{E}\).

Gamma = Matrix([A1, A2, A3]) * exp(I*kappa*x - I*omega*t)
E = Matrix([B1, B2, B3]) * exp(I*kappa*x - I*omega*t)

And compute the equations using the vector operators. Namely, the Laplace (vector.lap_vec()) and curl (vector.curl()) operators.

eq1 = c2**2 * lap_vec(Gamma) + K**2*curl(E) - Gamma.diff(t, 2)
eq2 = Q**2 * curl(Gamma) + c4**2*lap_vec(E) - 2*Q**2*E - E.diff(t, 2)
eq1 = simplify(eq1/exp(I*kappa*x - I*omega*t))
eq2 = simplify(eq2/exp(I*kappa*x  - I*omega*t))
eq = eq1.col_join(eq2)

We can compute the matrix for the system using .jacobian()

M = eq.jacobian([A1, A2, A3, B1, B2, B3])
M

\[\begin{split}\left[\begin{matrix}- c_{2}^{2} \kappa^{2} + \omega^{2} & 0 & 0 & 0 & 0 & 0\\0 & - c_{2}^{2} \kappa^{2} + \omega^{2} & 0 & 0 & 0 & - i K^{2} \kappa\\0 & 0 & - c_{2}^{2} \kappa^{2} + \omega^{2} & 0 & i K^{2} \kappa & 0\\0 & 0 & 0 & - 2 Q^{2} - c_{4}^{2} \kappa^{2} + \omega^{2} & 0 & 0\\0 & 0 & - i Q^{2} \kappa & 0 & - 2 Q^{2} - c_{4}^{2} \kappa^{2} + \omega^{2} & 0\\0 & i Q^{2} \kappa & 0 & 0 & 0 & - 2 Q^{2} - c_{4}^{2} \kappa^{2} + \omega^{2}\end{matrix}\right]\end{split}\]

And, we are interested in the determinant of the matrix \(M\).

factor(M.det())

\[\left(c_{2} \kappa - \omega\right) \left(c_{2} \kappa + \omega\right) \left(2 Q^{2} + c_{4}^{2} \kappa^{2} - \omega^{2}\right) \left(- K^{2} Q^{2} \kappa^{2} + 2 Q^{2} c_{2}^{2} \kappa^{2} - 2 Q^{2} \omega^{2} + c_{2}^{2} c_{4}^{2} \kappa^{4} - c_{2}^{2} \kappa^{2} \omega^{2} - c_{4}^{2} \kappa^{2} \omega^{2} + \omega^{4}\right)^{2}\]

The roots for this polynomial (in \(\omega^2\)) represent the dispersion relations.

disps = solve(M.det(), omega**2)
for disp in disps:
    display(disp)

\[c_{2}^{2} \kappa^{2}\]

\[2 Q^{2} + c_{4}^{2} \kappa^{2}\]

\[Q^{2} + \frac{c_{2}^{2} \kappa^{2}}{2} + \frac{c_{4}^{2} \kappa^{2}}{2} - \frac{1}{2} \sqrt{4 K^{2} Q^{2} \kappa^{2} + 4 Q^{4} - 4 Q^{2} c_{2}^{2} \kappa^{2} + 4 Q^{2} c_{4}^{2} \kappa^{2} + c_{2}^{4} \kappa^{4} - 2 c_{2}^{2} c_{4}^{2} \kappa^{4} + c_{4}^{4} \kappa^{4}}\]

\[Q^{2} + \frac{c_{2}^{2} \kappa^{2}}{2} + \frac{c_{4}^{2} \kappa^{2}}{2} + \frac{1}{2} \sqrt{4 K^{2} Q^{2} \kappa^{2} + 4 Q^{4} - 4 Q^{2} c_{2}^{2} \kappa^{2} + 4 Q^{2} c_{4}^{2} \kappa^{2} + c_{2}^{4} \kappa^{4} - 2 c_{2}^{2} c_{4}^{2} \kappa^{4} + c_{4}^{4} \kappa^{4}}\]

References

[NOWACKI]

Nowacki, W. (1986). Theory of asymmetric elasticity. pergamon Press, Headington Hill Hall, Oxford OX 3 0 BW, UK, 1986.

[GUARIN]

Guarín-Zapata, N., Gomez, J., Valencia, C., Dargush, G. F., & Hadjesfandiari, A. R. (2020). Finite element modeling of micropolar-based phononic crystals. Wave Motion, 92, 102406.

Continuum Mechanics package

Vector calculus module

vector.antisym_grad(A, coords=(x, y, z), h_vec=(1, 1, 1))

Antisymmetric part of the gradient of a vector function A.

A : Matrix (3, 1), list

Vector function to compute the gradient from.

coords : Tuple (3), optional

Coordinates for the new reference system. This is an optional parameter it takes (x, y, z) as default.

h_vec : Tuple (3), optional

Scale coefficients for the new coordinate system. It takes (1, 1, 1), as default.

antisym_grad: Matrix (3, 3)

Matrix with the components of the antisymmetric part of the gradient. The position (i, j) has as components diff(A[i], coords[j].

vector.biharmonic(u, coords=(x, y, z), h_vec=(1, 1, 1))

Bilaplacian of the scalar function u.

u : SymPy expression

Scalar function to compute the bilaplacian from.

coords : Tuple (3), optional

Coordinates for the new reference system. This is an optional parameters, and it takes a cartesian (x, y, z), as default.

h_vec : Tuple (3), optional

Scale coefficients for the new coordinate system. It takes (1, 1, 1), as default.

bilaplacian: Sympy expression

Bilaplacian of u.

vector.curl(A, coords=(x, y, z), h_vec=(1, 1, 1))

Curl of a function vector A.

A : Matrix, List

Vector function to compute the curl from.

coords : Tuple (3), optional

Coordinates for the new reference system. This is an optional parameter it takes (x, y, z) as default.

h_vec : Tuple (3), optional

Scale coefficients for the new coordinate system. It takes (1, 1, 1), as default.

curl : Matrix (3, 1)

Column vector with the curl of A.

vector.div(A, coords=(x, y, z), h_vec=(1, 1, 1))

Divergence of the vector function A.

A : Matrix, list

Vector function to compute the divergence from.

coords : Tuple (3), optional

Coordinates for the new reference system. This is an optional parameter it takes (x, y, z) as default.

h_vec : Tuple (3), optional

Scale coefficients for the new coordinate system. It takes (1, 1, 1), as default.

divergence: SymPy expression

Divergence of A.

vector.div_tensor(tensor, coords=(x, y, z), h_vec=(1, 1, 1))

Divergence of a (second order) tensor

tensor : Matrix (3, 3)

Tensor function function to compute the divergence from.

coords : Tuple (3), optional

Coordinates for the new reference system. This is an optional parameter it takes (x, y, z) as default.

h_vec : Tuple (3), optional

Scale coefficients for the new coordinate system. It takes (1, 1, 1), as default.

divergence: Matrix

Divergence of tensor.

[RICHARDS]

Rowland Richards. Principles of Solids Mechanics. CRC Press, 2011.

vector.dual_tensor(vec)

Compute the dual tensor for an axial vector

In index notation, the dual is defined by

\[C_{ij} = \epsilon_{ijk} C_k\]

where \(\epsilon_{ijk}\) is the Levi-Civita symbol.

vec : Matrix (3)

Axial vector.

dual: Matrix (3, 3)

Second order matrix that is dual of vec.

[ARFKEN]

George Arfken, Hans J. Weber and Frank Harris. Mathematical methods for physicists, Elsevier, 2013.

vector.dual_vector(tensor)

Compute the dual (axial) vector for an anti-symmetric tensor

In index notation, the dual is defined by

\[C_{i} = \frac{1}{2}\epsilon_{ijk} C_{jk}\]

where \(\epsilon_{ijk}\) is the Levi-Civita symbol.

tensor : Matrix (3, 3)

Second order tensor.

dual: Matrix (3)

Axial vector that is the dual of tensor.

[ARFKEN]

George Arfken, Hans J. Weber and Frank Harris. Mathematical methods for physicists, Elsevier, 2013.

vector.grad(u, coords=(x, y, z), h_vec=(1, 1, 1))

Compute the gradient of a scalara function phi.

u : SymPy expression

Scalar function to compute the gradient from.

coords : Tuple (3), optional

Coordinates for the new reference system. This is an optional parameters, and it takes a cartesian (x, y, z), as default.

h_vec : Tuple (3), optional

Scale coefficients for the new coordinate system. It takes (1, 1, 1), as default.

gradient: Matrix (3, 1)

Column vector with the components of the gradient.

vector.grad_vec(A, coords=(x, y, z), h_vec=(1, 1, 1))

Gradient of a vector function A.

A : Matrix (3, 1), list

Vector function to compute the gradient from.

coords : Tuple (3), optional

Coordinates for the new reference system. This is an optional parameter it takes (x, y, z) as default.

h_vec : Tuple (3), optional

Scale coefficients for the new coordinate system. It takes (1, 1, 1), as default.

gradient: Matrix (3, 3)

Matrix with the components of the gradient. The position (i, j) has as components diff(A[i], coords[j].

vector.lap(u, coords=(x, y, z), h_vec=(1, 1, 1))

Laplacian of the scalar function u.

u : SymPy expression

Scalar function to compute the laplacian from.

coords : Tuple (3), optional

Coordinates for the new reference system. This is an optional parameters, and it takes a cartesian (x, y, z), as default.

h_vec : Tuple (3), optional

Scale coefficients for the new coordinate system. It takes (1, 1, 1), as default.

laplacian: Sympy expression

Laplacian of u.

vector.lap_vec(A, coords=(x, y, z), h_vec=(1, 1, 1))

Laplacian of a vector function A.

A : Matrix, List

Vector function to compute the laplacian from.

coords : Tuple (3), optional

Coordinates for the new reference system. This is an optional parameter it takes (x, y, z) as default.

h_vec : Tuple (3), optional

Scale coefficients for the new coordinate system. It takes (1, 1, 1), as default.

laplacian : Matrix (3, 1)

Column vector with the components of the Laplacian.

vector.levi_civita(i, j, k)

Levi-Civita symbol

vector.scale_coeff(r_vec, coords)

Compute scale coefficients for the vector tranform given by r_vec.

r_vec : Matrix (3, 1)

Transform vector (x, y, z) as a function of coordinates u1, u2, u3.

coords : Tuple (3)

Coordinates for the new reference system.

h_vec : Tuple (3)

Scale coefficients.

vector.scale_coeff_coords(coord_sys, coords, a=1, b=1, c=1)

Return scale factors for predefined coordinate system.

coord_sys : string

Coordinate system.

coords : Tuple (3)

Coordinates for the new reference system.

a : SymPy expression, optional

Additional parameter for some coordinate systems.

b : SymPy expression, optional

Additional parameter for some coordinate systems.

c : SymPy expression, optional

Additional parameter for some coordinate systems.

h_vec : Tuple (3)

Scale coefficients.

[ORTHO]

Wikipedia contributors, ‘Orthogonal coordinates’, Wikipedia, The Free Encyclopedia, 2019

vector.sym_grad(A, coords=(x, y, z), h_vec=(1, 1, 1))

Symmetric part of the gradient of a vector function A.

A : Matrix (3, 1), list

Vector function to compute the gradient from.

coords : Tuple (3), optional

Coordinates for the new reference system. This is an optional parameter it takes (x, y, z) as default.

h_vec : Tuple (3), optional

Scale coefficients for the new coordinate system. It takes (1, 1, 1), as default.

sym_grad: Matrix (3, 3)

Matrix with the components of the symmetric part of the gradient. The position (i, j) has as components diff(A[i], coords[j].

vector.transform_coords(coord_sys, coords, a=1, b=1, c=1)

Return transformation for predefined coordinate systems.

coord_sys : string

Coordinate system.

coords : Tuple (3)

Coordinates for the new reference system.

a : SymPy expression, optional

Additional parameter for some coordinate systems.

b : SymPy expression, optional

Additional parameter for some coordinate systems.

c : SymPy expression, optional

Additional parameter for some coordinate systems.

h_vec : Tuple (3)

Scale coefficients.

[ORTHO]

Wikipedia contributors, ‘Orthogonal coordinates’, Wikipedia, The Free Encyclopedia, 2019

vector.unit_vec_deriv(vec_i, coord_j, coords=(x, y, z), h_vec=(1, 1, 1))

Compute the derivatives of unit vectors with respect to coordinates

The derivative is defined as

\[\begin{split}\frac{\partial{\hat{\mathbf{e}}_i}}{\partial u_j} = \begin{cases} \hat{\mathbf{e}_j} \frac{1}{h_i} \frac{\partial{h_j}}{\partial u_i} &\text{if } i\neq j\\ -\sum_{\substack{k=1\\ k\neq i}}^3 \hat{\mathbf{e}_k} \frac{1}{h_k} \frac{\partial{h_i}}{\partial u_k} &\text{if } i = j\\ \end{cases}\, ,\end{split}\]

as presented in [ARFKEN].

vec_i : int

Number of the unit vector (0, 1, 2).

coord_j : int

Number of the coordinate (0, 1, 2).

coords : Tuple (3), optional

Coordinates for the new reference system. This is an optional parameter it takes (x, y, z) as default.

h_vec : Tuple (3), optional

Scale coefficients for the new coordinate system. It takes (1, 1, 1), as default.

deriv: Matrix (3, 1)

Derivative of the i-th unit vector with respect to the j-th coordinate.

[ARFKEN]

George Arfken, Hans J. Weber and Frank Harris. Mathematical methods for physicists, Elsevier, 2013.

Solid mechanics module

solids.c_cst(u, parameters, coords=(x, y, z), h_vec=(1, 1, 1))

Corrected-Couple-Stress-Theory (C-CST) elasticity operator of a vector function u, as presented in [CST].

u : Matrix (3, 1), list

Vector function to apply the Navier-Cauchy operator.

parameters : tuple

Material parameters in the following order:

lamda : float

Lamé’s first parameter.

mu : float, > 0

Lamé’s second parameter.

eta : float, > 0

Couple stress modulus in C-CST.

coords : Tuple (3), optional

Coordinates for the new reference system. This is an optional parameter it takes (x, y, z) as default.

h_vec : Tuple (3), optional

Scale coefficients for the new coordinate system. It takes (1, 1, 1), as default.

c_cst : Matrix (3, 1)

Components of the C-CST operator applied to the displacement vector

[CST]

(1, 2, 3) Ali R. Hadhesfandiari, Gary F. Dargush. Couple stress theory for solids. International Journal for Solids and Structures, 2011, 48, 2496-2510.

solids.disp_def_cst(u, coords=(x, y, z), h_vec=(1, 1, 1))

Compute strain measures for C-CST elasticity, as defined in [CST].

u : Matrix (3, 1), list

Displacement vector function to apply the micropolar operator.

coords : Tuple (3), optional

Coordinates for the new reference system. This is an optional parameter it takes (x, y, z) as default.

h_vec : Tuple (3), optional

Scale coefficients for the new coordinate system. It takes (1, 1, 1), as default.

strain : Matrix (3, 3)

Strain tensor.

curvature : Matrix (3, 3)

Curvature tensor.

solids.disp_def_micropolar(u, phi, coords=(x, y, z), h_vec=(1, 1, 1))

Compute strain measures for micropolar elasticity, as defined in [NOW].

u : Matrix (3, 1), list

Displacement vector function to apply the micropolar operator.

phi : Matrix (3, 1), list

Microrrotation (pseudo)-vector function to apply the micropolar operator.

coords : Tuple (3), optional

Coordinates for the new reference system. This is an optional parameter it takes (x, y, z) as default.

h_vec : Tuple (3), optional

Scale coefficients for the new coordinate system. It takes (1, 1, 1), as default.

strain : Matrix (3, 3)

Strain tensor.

curvature : Matrix (3, 3)

Curvature tensor.

solids.micropolar(u, phi, parameters, coords=(x, y, z), h_vec=(1, 1, 1))

Micropolar operator of a vector function u, as defined in [NOW].

u : Matrix (3, 1), list

Displacement vector function to apply the micropolar operator.

phi : Matrix (3, 1), list

Microrrotation (pseudo)-vector function to apply the micropolar operator.

parameters : tuple

Material parameters in the following order:

lamda : float

Lamé’s first parameter.

mu : float, > 0

Lamé’s second parameter.

alpha : float, > 0

Micropolar parameter.

beta : float

Micropolar parameter.

gamma : float, > 0

Micropolar parameter.

epsilon : float, > 0

Micropolar parameter.

coords : Tuple (3), optional

Coordinates for the new reference system. This is an optional parameter it takes (x, y, z) as default.

h_vec : Tuple (3), optional

Scale coefficients for the new coordinate system. It takes (1, 1, 1), as default.

disp_op : Matrix (3, 1)

Displacement components.

rot_op : Matrix (3, 1)

Microrrotational components.

[NOW]

(1, 2) Witold Nowacki. Theory of micropolar elasticity. International centre for mechanical sciences, Courses and lectures, No. 25. Berlin: Springer, 1972.

solids.navier_cauchy(u, parameters, coords=(x, y, z), h_vec=(1, 1, 1))

Navier-Cauchy operator of a vector function u.

u : Matrix (3, 1), list

Vector function to apply the Navier-Cauchy operator.

parameters : tuple

Material parameters in the following order:

lamda : float

Lamé’s first parameter.

mu : float, > 0

Lamé’s second parameter.

coords : Tuple (3), optional

Coordinates for the new reference system. This is an optional parameter it takes (x, y, z) as default.

h_vec : Tuple (3), optional

Scale coefficients for the new coordinate system. It takes (1, 1, 1), as default.

navier_op : Matrix (3, 1)

Components of the Navier-Cauchy operator applied to the displacement vector.

solids.strain_stress(strain, parameters)

Return the stress tensor for a given strain tensor and material properties lambda and mu.

u : Matrix (3, 3)

Strain tensor.

parameters : tuple

Material parameters in the following order:

lamda : float

Lamé’s first parameter.

mu : float, > 0

Lamé’s second parameter.

disp_op : Matrix (3, 1)

Displacement components.

solids.strain_stress_cst(strain, curvature, parameters)

Return the force-stress and couple-stress tensor for given strain and curvature tensors and material parameters for the Corrected-Couple-Stress-Theory (C-CST), as presented in [CST].

strain : Matrix (3, 3)

Strain tensor.

curvature : Matrix (3, 3)

Curvature tensor.

parameters : tuple

Material parameters in the following order:

lamda : float

Lamé’s first parameter.

mu : float, > 0

Lamé’s second parameter.

eta : float, > 0

Couple stress modulus in C-CST.

solids.strain_stress_micropolar(strain, curvature, parameters)

Return the force-stress and couple-stress tensor for given strain and curvature tensors and material parameters.

strain : Matrix (3, 3)

Strain tensor.

curvature : Matrix (3, 3)

Curvature tensor.

parameters : tuple

Material parameters in the following order:

lamda : float

Lamé’s first parameter.

mu : float, > 0

Lamé’s second parameter.

alpha : float, > 0

Micropolar parameter.

beta : float

Micropolar parameter.

gamma : float, > 0

Micropolar parameter.

epsilon : float, > 0

Micropolar parameter.

Visualization module

Functions to aid the visualization of mathematical entities such as second rank tensors.

visualization.christofel_eig(C, nphi=30, nth=30)

Compute surfaces of eigenvalues for the Christoffel stiffness

For every direction

\[\mathbf{n} =(\sin(\theta)\cos(\phi),\sin(\theta)\sin(\phi),\cos(\theta))\]

computes the eigenvalues of the Christoffel stiffness tensor. These eigencalues are \(\rho v_p^2\), where \(\rho\) is the mass density and \(v_p\) is the phase speed.

C : (6,6) array

Stiffness tensor in Voigt notation.

nphi : (optional) int

Number of partitions in the azimut angle (phi).

nth : (optional) int

Number of partitions in the cenit angle (theta).

V1, V2, V3 : (nphi, nth) arrays

Eigenvalues for the desired discretization.

phi_vec : (nphi) array

Array with azimut angles.

theta_vec : (nth) array

Array with cenit angles.

visualization.mohr2d(stress, ax=None)

Plot Mohr circle for a 2D tensor

stress : ndarray

Stress tensor.

ax : Matplotlib axes, optional

Axes where the plot is going to be added.

[BRAN]

Brannon, R. (2003). Mohr’s Circle and more circles. Poslední revize, 29(10).

visualization.mohr3d(stress, ax=None)

Plot 3D Mohr circles

stress : ndarray

Stress tensor.

ax : Matplotlib axes, optional

Axes where the plot is going to be added.

visualization.plot_surf(R, phi, theta, title='Wave speed', **kwargs)

Plot the surface given by R(phi, theta).

R : (m,n) ndarray

Radius function.

phi_vec : (m,n) ndarray

Meshgrid for the azimut angle (phi).

theta_vec : (m,n) ndarray

Meshgrid for the cenit angle (theta).

**kwargs : keyword arguments (optional)

Keyword arguments for mlab.mesh.

surf : mayavi mesh

Mayavi mesh for the surface R(phi, theta).

visualization.traction_circle(stress, npts=48, ax=None)

Visualize a second order tensor as a collection of tractions vectors over a circle.

stress : ndarray

Stress tensor.

npts : int, optional

Number of vector to plot over the circle.

ax : Matplotlib axes, optional

Axes where the plot is going to be added.

Contributing

Contributions are welcome, and they are greatly appreciated! Every little bit helps, and credit will always be given.

You can contribute in many ways:

Types of Contributions

Report Bugs

Report bugs at https://github.com/nicoguaro/continuum_mechanics/issues.

If you are reporting a bug, please include:

• Your operating system name and version.
• Any details about your local setup that might be helpful in troubleshooting.
• Detailed steps to reproduce the bug.

Fix Bugs

Look through the GitHub issues for bugs. Anything tagged with “bug” and “help wanted” is open to whoever wants to implement it.

Implement Features

Look through the GitHub issues for features. Anything tagged with “enhancement” and “help wanted” is open to whoever wants to implement it.

Write Documentation

Continuum Mechanics could always use more documentation, whether as part of the official Continuum Mechanics docs, in docstrings, or even on the web in blog posts, articles, and such.

Submit Feedback

The best way to send feedback is to file an issue at https://github.com/nicoguaro/continuum_mechanics/issues.

If you are proposing a feature:

• Explain in detail how it would work.
• Keep the scope as narrow as possible, to make it easier to implement.
• Remember that this is a volunteer-driven project, and that contributions are welcome :)

Get Started!

Ready to contribute? Here’s how to set up continuum_mechanics for local development.

1. Fork the continuum_mechanics repo on GitHub.

2. Clone your fork locally:

   $ git clone git@github.com:your_name_here/continuum_mechanics.git
   

3. Install your local copy into a virtualenv. Assuming you have virtualenvwrapper installed, this is how you set up your fork for local development:

   $ mkvirtualenv continuum_mechanics
   $ cd continuum_mechanics/
   $ python setup.py develop
   

4. Create a branch for local development:

   $ git checkout -b name-of-your-bugfix-or-feature
   

   Now you can make your changes locally.

5. When you’re done making changes, check that your changes pass flake8 and the tests, including testing other Python versions with tox:

   $ flake8 continuum_mechanics tests
   $ python setup.py test or py.test
   $ tox
   

   To get flake8 and tox, just pip install them into your virtualenv.

6. Commit your changes and push your branch to GitHub:

   $ git add .
   $ git commit -m "Your detailed description of your changes."
   $ git push origin name-of-your-bugfix-or-feature
   

7. Submit a pull request through the GitHub website.

Pull Request Guidelines

Before you submit a pull request, check that it meets these guidelines:

1. The pull request should include tests.
2. If the pull request adds functionality, the docs should be updated. Put your new functionality into a function with a docstring, and add the feature to the list in README.rst.
3. The pull request should work for Python 3.7 to 3.9, and for PyPy. Check https://travis-ci.org/nicoguaro/continuum_mechanics/pull_requests and make sure that the tests pass for all supported Python versions.

Tips

To run a subset of tests:

$ py.test tests.test_continuum_mechanics

Deploying

A reminder for the maintainers on how to deploy. Make sure all your changes are committed (including an entry in HISTORY.rst). Then run:

$ bumpversion patch # possible: major / minor / patch
$ git push
$ git push --tags

Travis will then deploy to PyPI if tests pass.

Credits

Development Lead

• Nicolás Guarín-Zapata @nicoguaro

Contributors

None yet. Why not be the first?

History

0.1.0 (2018-04-21)

• First release on PyPI.

0.1.2 (2019-04-04)

• Add elasticity operators.

0.2.0 (2019-12-02)

• Python 2 support dropped.
• Fix vector gradient for curvilinear coordinates.
• Add divergence for second-order tensors.

0.2.2 (2021-04-12)

• Python 3.6 support dropped.

Indices and tables

• Index
• Module Index
• Search Page