nicoguaro / continuum_mechanics Goto Github PK

The gradient of a vector returns wrong results for coordinates different from Cartesian.

• Continuum Mechanics version: 0.3.0

Description

The rotation matrix given in http://solidmechanics.org/Text/Chapter3_2/Chapter3_2.php under section 3.2.11 is not the same as what I get from the package. It might be an issue with clockwise and counterclockwise rotation, as the numbers are the same, the signs are not always.

As far as I understand it, a counterclockwise rotation should be positive about the rotation axis according to the right-hand rule. This is also what's shown in his book and the figure under 3.2.11 on the webpage.

I know this is a version under development, but I checked your code and did not find any error.

What I Did

Computed the two rotation matrices as below. They don't match.

theta = sym.pi/6
c = sym.cos(theta)
s = sym.sin(theta)

# counterclockwise, z-axis, 30°
# Taken from https://en.wikipedia.org/wiki/Rotation_matrix#Basic_rotations
Q = sym.Matrix([[c, -s, 0], [s, c, 0], [0, 0, 1]])     

M = sym.zeros(6)
for i in range(3):
    for j in range(3):
        M[i, j] = Q[i, j]**2
        M[i, j + 3] = 2 * Q[i, (j + 1) % 3] * Q[i, (j + 2) % 3]
        M[i + 3, j] = Q[(i + 1) % 3, j] * Q[(i + 2) % 3, j]
        M[i + 3, j + 3] = Q[(i + 1) % 3, (j + 1) % 3] * \
            Q[(i + 2) % 3, (j + 2) % 3] + \
            Q[(i + 1) % 3, (j + 2) % 3] * \
            Q[(i + 2) % 3, (j + 1) % 3]

# Counterclockwise, z-axis
M_ref = np.array([[c**2, s**2, 0, 0, 0, 2*c*s],
                  [s**2, c**2, 0, 0, 0, -2*c*s],
                  [0, 0, 1, 0, 0, 0],
                  [0, 0, 0, c, s, 0],
                  [0, 0, 0, -s, c, 0],
                  [-c*s, c*s, 0, 0, 0, c**2-s**2]])

print(M_ref == np.array(M))

Output:

[[ True  True  True  True  True False]
 [ True  True  True  True  True False]
 [ True  True  True  True  True  True]
 [ True  True  True  True  True  True]
 [ True  True  True  True  True  True]
 [False False  True  True  True  True]]

Depending on the context people use different combinations of the material constants in elasticity. We need a way to get any combination of two of them from any other two, as shown in:

https://en.wikipedia.org/wiki/Lam%C3%A9_parameters

Regarding Mohr circles in 2D it would be useful to have optional parameters to turn off some of the details in the graphics, such as the center or the line connecting the current state.

Add coordinate transformations for the supported coordinate systems.

• Continuum Mechanics version: 0.2.2
• Python version: 3.7.10

Description

The divergence of the stress tensor is not 0 while the Navier-Cauchy operator is zero for the displacement vector.

What I Did

Paste the command(s) you ran and the output.
If there was a crash, please include the traceback here.

from sympy import *
from continuum_mechanics.vector import div_tensor
from continuum_mechanics.solids import navier_cauchy

r, theta, phi = symbols("r theta phi")
T, mu = symbols("T mu")

h = (1, r, r*sin(theta))
coords = (r, theta, phi)
disp = Matrix([[0], [0], [T*sin(theta)/(8*pi*mu*r**2)]])
stress = Matrix([[0, 0, -3*T*sin(theta)/(8*pi*r**3)], [0, 0, 0], [-3*T*sin(theta)/(8*pi*r**3), 0, 0]])

And

navier_cauchy(disp, [1, mu], coords, h)

gives Matrix([[0], [0], [0]]), while

div_tensor(stress, coords, h)

gives Matrix([[0], [0], [3*T*sin(theta)/(8*pi*r**4)]]).

We can add Airy stress functions to the package since the vector operators are implemented. A good starting point is given in the following notebook:

https://github.com/nicoguaro/notebooks_examples/blob/master/Airy%20stress%20functions.ipynb

It could be possible to add Michelle general solution for 2D elasticity in polar coordinates.

• J. H. Michell, 1899, On the direct determination of stress in an elastic solid,
  with application to the theory of plates, Proceedings of the London Mathematical
  Society, vol. 31, pages 100-124.