Make the 2D formulations consistent with plane strain assumption

[YiminJin]

I found that the 2D formulations in ASPECT are not fully consistent with the plane strain assumption when trying to add plastic dilation to the Drucker-Prager rheology (#6373 ). The inconsistency is caused by two functions frequently used by the material models:

1. deviator(). This function is defined by $\text{dev}(A) = A - \dfrac{1}{\text{dim}}\text{trace}(A)I$, while under the plane strain assumption the deviatoric strain rate should be $\varepsilon' = \varepsilon - \dfrac{1}{3}\text{trace}(\varepsilon)I$.
2. second_invariant(). The second invariant of a deviatoric symmetric tensor can be expressed as $II_A = -\dfrac{1}{2}A_{ij}A_{ij}$. However, under the plane strain assumption, the second invariant of deviatoric strain rate should be $II_\varepsilon = -\dfrac{1}{2}[\varepsilon'^2_{11} + \varepsilon'^2_{22} + \varepsilon'^2_{33} + 2\varepsilon'^2_{12}]$.

I have tried the plane strain version of the two functions in the strip footing test and Kaus' [2010] test with associated plastic flow, and a big improvement on the sharpness of the shear bands can be observed from the results (left: before the modification; right: after the modification):

However, I am not sure if deviator() and second_invariant() are the only functions inconsistent with the plane strain assumption. Does anyone see something that I have missed?

[bangerth]

Is the difference due to the question of whether we assume plane strain vs plane stress? What does the manual say we use?

[YiminJin]

Is the difference due to the question of whether we assume plane strain vs plane stress? What does the manual say we use?

I think the manual does not mention the term "plane strain". There is a subsection that explains why we use 1/3 in the system matrix in 2D (https://aspect-documentation.readthedocs.io/en/latest/user/methods/basic-equations/2d-models.html), which consists with plane strain assumption.

[bobmyhill]

@YiminJin thanks for opening this Issue!

@bangerth This looks like a bug that we should fix, but we should definitely all agree before we make a decision, because it is a fairly major breaking change. It's fairly clear in the manual that we intend to use plane strain in 2D models. In incompressible models, the distinction is unimportant. For example, $\text{trace}(\varepsilon) = 0$ in 2D or 3D-plane strain, and so the expression for the deviator is insensitive to whether a factor 1/2 or 1/3 is used. In compressible models undergoing plane strain, $du_y / dy \neq -du_x / dx$, and so $\text{trace}(\varepsilon) \neq 0$ is not equal to zero.

I think the plane strain approximation should also produce a plane stress. I think an important point is that in compressible plane strain, $\sigma_{zz}$ should change rapidly where the material is being compressed or dilated, and this can drive localised plastic failure.

Do we have any 2.5D compressible tests where the model is 1 cell wide? I think this would be a good test of behaviour.

[YiminJin]

Do we have any 2.5D compressible tests where the model is 1 cell wide? I think this would be a good test of behaviour.

I am going to try that out.

After discussion with @naliboff and @gassmoeller , we decided to limit #6373 to modifications on plastic dilation only, and leave the plane strain issues to a following-up PR. So I think it would be OK to test the 2.5D case with #6373 and see what the results look like.

[bangerth]

@bobmyhill I had not thought about the fact that if a material is incompressible, plane strain = plane stress. That's a good point. I think we ought to revisit the discussion of what 2d actually means in the manual. I might put that on my list of things to do this afternoon. I do recall, too, that we intended to interpret 2d as "plane strain".

[YiminJin]

I think the plane strain approximation should also produce a plane stress.

I think so too. We considered about this idea in the discussion, and we thought it is possibly right because it explains why the inconsistency has not led to invalid results before (there are few benchmarks with both compressibility and strain concentration).

[bangerth]

See #6459 for some additional comments.

[bobmyhill]

Thanks @bangerth , from #6459 I realise that I misunderstood the meaning of plane stress (not just $\sigma_{xz}=\sigma_{yz}=0$ and $\sigma_{zz}=c$, but rather $\sigma_{xz}=\sigma_{yz}=\sigma_{zz}=0$).

@YiminJin that's a good point. Also, I think the anisotropic strain due to compression and heating is (in most places where it matters) slow compared with convective strain and viscous relaxation, moderating the effect of applying a different factor to the trace of the strain rate.

[YiminJin]

I carried out experiments with the strip-footing test (with one level coarser than the figures posted in the first comment). Here are the results of the strain rate field:

The 2.5D results are very similar to the 2D results under the plane strain assumption, which implies that the modifications on deviator() and second_invariant() works well for viscoelastic plastic problems. I am going to pull a request for those modifications.

[bobmyhill]

This is brilliant, thank you for working on this @YiminJin!

When you open the PR, it would be great to have a really simple test to show the difference between the old and new models. Perhaps a simple downwards translation of compressible material under gravity? I guess the stress evolution will differ between the current code and your modifications.

[YiminJin]

@bobmyhill Thank you for the suggestion. Could you describe the model with more details? Or is there a prm file for a similar model? (I do not quite understand the meaning of "downwards translation").

[bobmyhill]

The model I had in mind is a 1.5D box, lets say 1000 km deep x 1 cell wide. Gravity is Earth-like. The sides are free slip, the top has prescribed inflow and the bottom allows outflow. The material inside is compressible, with an isentropic temperature gradient, and no shear heating.

As material moves downwards, it should compress under gravity, and that compression under plane-strain should generate a deviatoric stress. A balance will be reached between the stress generated by downward motion and viscous relaxation.

Does that make sense? I'm sure there are loads of other simple tests that would also highlight the discrepancies.