Problem with maximum horizontal compressive stress computation

I have been looking at the math behind the maximum horizontal compression post-processor, because I would like to change it to a maximum horizontal shear stress post-processor. Looked at the code and the manual, and I was left with a factor two in de derivation. At this point the equations are still the same:

f(\alpha) = \mathbf n^T \sigma_c \mathbf n = (\mathbf u^T \sigma_c \mathbf u)(\cos\alpha)^2 + 2(\mathbf u^T \sigma_c \mathbf v)(\cos\alpha)(\sin\alpha) +(\mathbf v^T \sigma_c \mathbf v)(\sin\alpha)^2

I do a substitution to simplify:

a = \mathbf u^T \sigma_c \mathbf u \\ b = \mathbf u^T \sigma_c \mathbf v \\ c = \mathbf v^T \sigma_c \mathbf v \\ f(\alpha) = a\cos^2\alpha + 2b(\cos\alpha)(\sin\alpha) +c\sin^2\alpha

The derivative of this is (see for example wolfram alpha):

f'(\alpha) = a*-2\sin(\alpha)\cos(\alpha)+2b(\cos^2(\alpha)-\sin^2(\alpha))+c*2\sin(\alpha)\cos(\alpha)

I will use the following trigonomic identities from wikipedia:

\tan(2\alpha) = \frac{2\tan(\alpha)}{1-\tan^2(\alpha)} \\ \sin(2\alpha) = 2\sin(\alpha)\cos(\alpha) = \frac{2\tan(\alpha)}{1+\tan^2(\alpha)} \\ \cos(2\alpha) = \cos^2(\alpha)-\sin^2(\alpha) = \frac{1-\tan^2(\alpha)}{1+\tan^2(\alpha)} \\

We use these idententies to substitute:

f'(\alpha) = -a*2\sin(\alpha)\cos(\alpha)+2b(\cos^2(\alpha)-\sin^2(\alpha))+c*2\sin(\alpha)\cos(\alpha) \\ = -a*\frac{2\tan(\alpha)}{1+\tan^2(\alpha)}+2b(\frac{1-\tan^2(\alpha)}{1+\tan^2(\alpha)})+c*\frac{2\tan(\alpha)}{1+\tan^2(\alpha)} \\ = \frac{1}{1+\tan^2(\alpha)}\left(-a*2\tan(\alpha)+2b\left(1-\tan^2(\alpha)\right)+c*2\tan(\alpha)\right) \\ = \frac{1-\tan^2(\alpha)}{1+\tan^2(\alpha)}\left(-a\frac{2\tan(\alpha)}{1-\tan^2(\alpha)}+2b+c \frac{2\tan(\alpha)}{1-\tan^2(\alpha)}\right) \\ = \frac{1-\tan^2(\alpha)}{1+\tan^2(\alpha)}\left(-a\tan(2\alpha)+2b+c \tan(2\alpha)\right) \\

Because we try to find alpha where the derivative is zero we say

-a\tan(2\alpha)+2b+c \tan(2\alpha) = 0 \\ a\tan(2\alpha)-c \tan(2\alpha) = 2b \\ \tan(2\alpha) = \frac{2b}{a-c} \\ \tan(2\alpha) = \frac{2\mathbf u^T \sigma_c \mathbf v}{\mathbf u^T \sigma_c \mathbf u-\mathbf v^T \sigma_c \mathbf v} \\

While the plugin description in the manual (which is also implemented as far as I can see) says:

\tan(2\alpha) = \frac{\mathbf u^T \sigma_c \mathbf v}{\mathbf u^T \sigma_c \mathbf u-\mathbf v^T \sigma_c \mathbf v} \\

Did I make a mistake in the derivation or is the description in the manual (and maybe the code) incorrect?

Hi Menno,
I would need to go through the derivation (I was not involved in writing this), and that will likely take a bit. Can you think of a simple testcase that would check which formula is correct? I.e. is there a case where the maximum compressive (or shear) stress should be in exactly one direction, and then check if the postprocessor computes it correctly?

Hey Rene,
Thanks for your reply. I have been thinking of doing this, because it would be good to have such test cases anyway, but I havn’t got around to thinking up a good case and testing it.

Hi Menno, Hi René,

Menno, your derivation looks correct to me. The only thing I’d say is that splitting the ellipsoid algebra into individual components obfuscates the derivation a bit. If you want to generalize to calculating normal and shear stresses on arbitrarily oriented planes, I think it would be simpler to stick with a pure tensor representation.

The current postprocessor should calculate the maximum (or minimum, depending on convention) eigenvector of the “horizontal stress ellipse”. This ellipse is the projection of the stress ellipsoid onto the horizontal plane. I’ll briefly rephrase the problem in mathematical terms:

If S is a symmetric m x m matrix representing an ellipsoid in an m-dimensional space, and C is an m x n orthonormal basis for an n-dimensional subspace, then the projection P of S onto that subspace is P = C^TSC. The eigenvectors and eigenvalues of P can be found in a number of ways (already implemented in deal.II/aspect).

The trigonometry in your derivation is the 2-dimensional variant of Jacobi’s method for finding eigenvalues of a matrix. Very briefly, the method finds the rotation angle \alpha which diagonalises the matrix [[a, b], [b, c]] by solving \tan(2\alpha) = 2b/(a-c), exactly the same as in your derivation.

I hope that helps :slight_smile:

Hey Bob and Rene,

First of all, thanks Bob for your reply and the generalization! It helped me quite a bit :slight_smile:

I have worked with Cedric (big thanks to him!) on a test example to see whether the 2 should or shouldn’t be there.

Let’s start with a 3d system in 2d simple shear:

\sigma = \begin{pmatrix}0 & \tilde\sigma & 0 \\ \tilde\sigma & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}

Since we can choose any u and v, let’s start with

u = e_x = \begin{pmatrix}1 \\ 0 \\ 0 \end{pmatrix} \text{ and } v = e_y = \begin{pmatrix}0 \\ 1 \\ 0 \end{pmatrix}

Given these vectors we expect that we need to rotate u and v 45 degrees, because of the state of stress. The tilde 2 indicates that the 2 might be there or not.

\tan(2\alpha) = \frac{\tilde 2 u^T \sigma v}{u^T \sigma u - v^T \sigma v} \\ = \frac{\tilde 2 \begin{pmatrix}1 \\ 0 \\ 0 \end{pmatrix}^T \begin{pmatrix}0 & \tilde \sigma & 0 \\ \tilde\sigma & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix}0 \\ 1 \\ 0 \end{pmatrix}}{\begin{pmatrix}1 \\ 0 \\ 0 \end{pmatrix}^T \begin{pmatrix}0 & \tilde\sigma & 0 \\ \tilde\sigma & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix}1 \\ 0 \\ 0 \end{pmatrix} - \begin{pmatrix}0 \\ 1 \\ 0 \end{pmatrix}^T \begin{pmatrix}0 & \tilde\sigma & 0 \\ \tilde\sigma & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}\begin{pmatrix}0 \\ 1 \\ 0 \end{pmatrix}} \\ = \frac{\tilde 2 \sigma }{0 - 0} \\ = \inf \\ \alpha = 0.5*\arctan(\inf) = \pi/4 = 45 \deg

Which is what we expect for this case, but doesn’t test whether the factor 2 is correct. Now let’s take two slightly different vectors:

u = \begin{pmatrix}\frac{1}{2} \\ 1 \\ 0 \end{pmatrix} \text{ and } v = \begin{pmatrix}-1 \\ \frac{1}{2} \\ 0 \end{pmatrix}

The angle between these vectors and the axes are (between v and e_x)

\theta = \arctan(\frac{1}{\frac{1}{2}}) = \arctan(2) \approx 63.4349 \deg \\ \alpha = 45-\theta = -18.4349 \deg

So we expect alpha to be 18.4349 degrees.

\tan(2\alpha) = \frac{\tilde 2 u^T \sigma v}{u^T \sigma u - v^T \sigma v} \\ = \frac{\tilde 2 \begin{pmatrix}\frac{1}{2} \\ 1 \\ 0 \end{pmatrix}^T \begin{pmatrix}0 & \tilde\sigma & 0 \\ \tilde\sigma & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix}-1 \\ \frac{1}{2} \\ 0 \end{pmatrix} }{\begin{pmatrix}\frac{1}{2} \\ 1 \\ 0 \end{pmatrix}^T \begin{pmatrix}0 & \tilde\sigma & 0 \\ \tilde\sigma & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix}\frac{1}{2} \\ 1 \\ 0 \end{pmatrix} - \begin{pmatrix}-1 \\ \frac{1}{2} \\ 0 \end{pmatrix} ^T \begin{pmatrix}0 & \tilde\sigma & 0 \\ \tilde\sigma & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix}-1 \\ \frac{1}{2} \\ 0 \end{pmatrix} } \\ = \frac{\tilde 2 \frac{3}{4}\sigma}{\sigma - \sigma} \\ = \frac{\tilde 2 \frac{3}{4}\sigma}{2 \sigma} \\ = \tilde 2 \frac{3}{8} \\ \alpha_1 = 0.5\arctan(-\frac{3}{8}) = -10,2750 \\ \alpha_2 = 0.5\arctan(-2\frac{3}{8}) = -18.4349

This shows that we need the 2 to be there.

Hi Menno, Hi Bob,
cool, this makes sense. Can you request the change and maybe add this as a testcase? I just checked and there does not seem to be a test for this postprocessor, I guess this again emphasizes to never trust untested code (also undermines our claim that everything in ASPECT is tested :frowning: ). Could you also add a changelog that this was changed?

Thanks for bringing this up :+1:.

Menno,
I’m to blame – I wrote this postprocessor in a morning while visiting Sarah a
couple of years ago. I must have made the mistake with the missing factor of 2
at the time. If you’re confident that it is really missing, then I have no
doubt that it really is.

The patch that added the original implementation is here:


I could have sworn that I added a testcase at the time, but apparently I
didn’t. But if you look at the picture in the PR, you can probably see how one
could construct one.

Best
W.

@gassmoeller I will make a pull request for the change. I just want to test a few more things.

@bangerth thanks for the link. I think I will use a slightly simpler setup of simple shear in x and y plus a z component, so that all directions are tested.