Reduce grain size due to plastic strain

[jdannberg]

In #5708, we added plasticity to the grain size model, but it does not yet affect grain size reduction.
This PR is a first attempt of making this happen. It assumes that as soon as material is yielding, all strain is plastic strain (which is what is assumed in other parts of ASPECT as well (and I think this is based on the idea that soon as material is yielding, we break the rock).

@bobmyhill Do you think this is reasonable? I started to think about how to partition the strain rate between diffusion, dislocation and plastic, had some ideas (which would all be very cmplicated and slow), discussed it with @gassmoeller and @MFraters, and this is what Menno suggested would be a good approximation.

@bobmyhill I also had no idea what would be a reasonable parameter for the Work fraction for boundary area change when yielding, do you happen to have an idea what could be reasonable? (Or is the answer just: we don't know?)

I also have not made any tests yet (and I'll do that as soon as we agree this is a reasonable way of implementing it).
For now, I just want to make sure that everything else still works the same way as before.

A part of this PR is also just refactoring things in the grain size material mode to avoid code duplication, and I am happy to make this a separate PR; let me know. I just through it would be good to also see the other changes so that it's clear why the refactoring is needed.

• I have followed the instructions for indenting my code.

For new features/models or changes of existing features:

• I have tested my new feature locally to ensure it is correct.
• I have created a testcase for the new feature/benchmark in the tests/ directory.
• I have added a changelog entry in the doc/modules/changes directory that will inform other users of my change.

[bobmyhill]

@bobmyhill Do you think this is reasonable?

• Short answer: Yes, I think that assuming that plasticity is the only mechanism that is active during brittle yielding is reasonable.
• Longer answer: There's some nuance in connecting what is happening in nature at the atom-scale with what we can reasonably simulate in a continuum model on m-km length scales.
  • Faulting involves discontinuities, so it is distinct from plasticity.
  • In our ASPECT we approximate brittle faulting as plastic. Soon we will add porosity generation into plasticity (via a dilation term), which will get closer to brittle behaviour, but it still won't mirror nature in some ways that matter for grain size reduction.
  • During faulting, new grain boundaries form by accumulation of dislocations, so in this sense fault-comminution is similar to grain size reduction by dislocation creep.
  • Stress fields and stress-strain relationships are very different under brittle failure compared with dislocation creep (e.g. shock loading), and so it makes sense to separate these in a modelling framework. See, for example: https://www.researchgate.net/profile/Charles-Sammis/publication/228515462_Mechanics_of_grain-size_reduction_in_fault_zones/links/00b49523af7c991e31000000/Mechanics-of-grain-size-reduction-in-fault-zones.pdf

@bobmyhill I also had no idea what would be a reasonable parameter for the Work fraction for boundary area change when yielding, do you happen to have an idea what could be reasonable? (Or is the answer just: we don't know?)

Short answer: Sorry, I don't know.
Longer answer: Last time I looked into this (10+ years ago) there were arguments about amount of energy that went into sensible and latent heat / damage / elastic or seismic energy / dilation (PV work). The newer literature will probably give us better estimates / arguments, but I haven't yet looked. There are a lot of field observations of grain size distributions in fault zones that we could try to match, but are we looking to match those observations?

I started to think about how to partition the strain rate between diffusion, dislocation and plastic, had some ideas (which would all be very cmplicated and slow)

I've also been thinking about this the last couple of days. The change in strain rate fraction between different viscous elements (solid lines in the figure below) assuming an isostress condition (i.e., what we do in diffusion_dislocation.cc) is beautifully modelled by a simple power-law function (dashed line, $\sum a / (1 + (\dot{\varepsilon}_0/\dot{\varepsilon}_T)^{1/c})^{d}) + f$, as found using a project I've adopted: https://github.com/equations-project/pyeq3):

but the integral of this expression (which would yield the strain rate fractions) involves hypergeometric functions that are horribly expensive to calculate.

I'll keep on thinking about this, but maybe the ``cheap'' way done by most other codes (harmonic viscosity averaging) is good enough?

[jdannberg]

@bobmyhill Thank you, that's super helpful!

For the work fraction for boundary area change when yielding, I just thought it would save me a lot of work if you had a good suggestion, but I'll have a look into the literature. For testing, I would use a constant value anyway and it does not really matter what that is. I don't think we need to match anything, the default should just be something that's somewhat reasonable. For now we could even just set it to 0.1, which is the same as the default for the work fraction for boundary area change for dislocation creep.

[bobmyhill]

For the work fraction for boundary area change when yielding, I just thought it would save me a lot of work

Yes, that's probably true. I had hoped that recent citations of Kanamori 2001 (https://www.google.co.uk/books/edition/Earthquake_Thermodynamics_and_Phase_Tran/z0x5Z4EjP6QC?hl=en&gbpv=1&dq=energy+seismic+efficiency&pg=PA293&printsec=frontcover) might help, but I didn't find much - understandably, most of the discussion revolves around the seismic efficiency, not the ratio of heat to other dissipative mechanisms (e.g. https://gfzpublic.gfz-potsdam.de/rest/items/item_108170/component/file_364681/content). We might need to ask someone closer to the field.

[gassmoeller]

Just some initial comments (I know this is WIP). Code structure looks good.

[gassmoeller]

wouldnt you need to check for if (stress > yield_stress) before these two lines? At the moment it looks like it is always yielding.

[gassmoeller]

Suggested change

	"The fraction $\\chi$ of work done by plastic tielding to change the grain boundary area. "
	"The fraction $\\chi$ of work done by plastic yielding to change the grain boundary area. "

[gassmoeller]

So when setting the new parameter to 0 nothing changes, which is good. But we also need a test that uses the new feature.