-
Notifications
You must be signed in to change notification settings - Fork 235
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
issue in morency_doin.cc when calculate the second invariant of strain rate #527
Comments
Dear Qingwen Zhang, I have no experience with the Morency Doin model, so I'll leave that Best wishes, @bobmyhill https://github.com/bobmyhill I recently notice something wrong — |
Yep, this still needs more testing, but it looks like the paper actually uses even a different definition for their second invariant (maybe incorrect, I'm not sure). The deal.ii second_invariant function calcuates something like: But what they call second invariant in the paper is defined as: This will be fixed soon. Sorry for the inconvenience. |
Hi again, Thanks for sending us your test results! Jonathan and I will fix these two Best wishes,
|
All, just as a point of reference: it surprises me that anyone would define the second invariant involving the square root. The first, second, and third invariants of a tensor are commonly defined so that they are linear, quadratic, and cubic functions of the tensor elements. See, for example, http://en.wikipedia.org/wiki/Cauchy_stress_tensor#Principal_stresses_and_stress_invariants and http://en.wikipedia.org/wiki/Invariants_of_tensors . If Morency and Doin really use a definition that involves the square root, can you please leave a comment in the code making it clear that their definition deviates from the common one, to avoid confusion by readers of the code at a later time? |
Yes, this definition and implementation per se is correct. I will recheck where the problem is as soon as possible. |
I had been sort of focused on this problem with e[0][0]*e[1][1] - e[0][1]^2 vs e[0][0]^2 + e[0][1]^2 but the square root thing is weirder I guess. @QwZhang do you know where this definition comes from? And yes, once this is all sorted out I will document this awkward definition in the comments and in the manual cookbook (which I am also working on improving). After this week I will finally have time to follow through with those fixes. |
For the record though, I am still kind of in favor of removing this material model from Aspect because I don't use it, and can't really vouch for its usefulness in most models. I think that the cookbook is instructive, but the material model should maybe only be distributed as part of the cookbook, and not one of the built-in material models because it is so specific to this one paper's model setup and I'm not convinced that it's appropriate for general purpose use. |
I think I find the problems.
So, for the following lines of implementation of plastic rheology: const SymmetricTensor<2,dim> strain_rate_dev = deviator(strain_rate);
double strain_rate_dev_inv = std::sqrt(0.5) * strain_rate_dev.norm();
strength = (std::max(pressure,0.0) * std::sin(phi) + cohesion * std::cos(phi));
const double viscosity = strength / (2.0 * strain_rate_dev_inv); the final result is correct but this strain_rate_dev_inv here actually equals to the square root of second invariant of deviatoric stress rate tensor, not the second invariant of deviatoric stress rate tensor as indicated by its name. Therefore, it's not very intuitive and can be replaced by: const SymmetricTensor<2,dim> strain_rate_dev = deviator(strain_rate);
const double strain_rate_dev_inv = second_invariant(strain_rate_dev);
strength = (std::max(pressure,0.0) * std::sin(phi) + cohesion * std::cos(phi));
const double viscosity = strength / (2.0 * std::sqrt(strain_rate_dev_inv)); This actually yield same results as shown in the figures above. For viscocity, the effective viscosity can be calculated in various way, in which both @jperryhouts, perhaps I have been misleading you, sorry for the inconvenience! Regards, |
On 06/11/2015 11:24 PM, Jonathan Perry-Houts wrote:
If div(u)=0, then e[0][0]=-e[1][1], so indeed the sign appears wrong.
Thanks. |
On 06/11/2015 11:30 PM, Jonathan Perry-Houts wrote:
I like the idea of having many material models. Nobody is forced to use them, |
@bobmyhill I recently notice something wrong in aspect/source/material_model/morency_doin.cc when it deals with the second invariant of strain rate: the final value of strain rate should be taken the square root if one have used deal.ii's function second_invariant(), in which this operation is not included internally. In short, it should be somthing like this:
not:
Or do I miss something that I don't knew?
The text was updated successfully, but these errors were encountered: