Implemented Associative viscoplastic flow rule with MCC model

[atpellet]

Implemented a simple associative viscoplastic flow rule for the MCC model.

In the original case of plasticity, the strain computed in src/simulation/plasticity.cpp, initially

$$\boldsymbol\varepsilon^{E, n+1} = \boldsymbol\varepsilon^{E, trial} - \dot\gamma \frac{\partial f}{\partial \boldsymbol\tau} \Delta t$$

becomes, with the associative viscosity

$$\boldsymbol\varepsilon^{E, n+1} = \boldsymbol\varepsilon^{E, trial} - \frac{\dot\gamma}{\eta} \frac{\partial f}{\partial \boldsymbol\tau} \Delta t$$

with $\eta$ the viscosity.

This modification has necessitated to create a new version of the MCCRMAExcplicit() in src/plasticity_helpers/mcc_rma_explicit.cpp in order the extract $\Delta\Gamma=\Delta t\dot\gamma$ with $\dot\gamma$ the plastic multiplier, such that $\boldsymbol{l}^P = \dot\gamma \frac{\partial y}{\partial \boldsymbol\tau}$.
In the code, $\Delta\Gamma$ is called Delta_GAMMA.

[larsblatny]

Hi @atpellet. This is fantastic, thanks for contributing your model!
Before merging, i just want to quickly review the code.

One first question: It seems that the new plasticity helper function only differs from MCCRMAExcplicit by making the delta_gamma returnable? Perhaps we can find a non-intrusive way to add this directly to the existing code so we don't need to duplicate the entire code for such a small change?

On another note, I have not seen this kind of viscoplasticity before so if you have some references that use or introduce similar models I would be very interested. One concern I have is in line 653 and 654 of plasticity.cpp where dydp and dydq (which you need to shift the (p,q)-point from the yield surface to a point away from the yield surafce given by eta) are calculated based on the point on the yield and not on the n+1 point. Are you sure this would give the same result? I suppose this could be easily verified. It is also clear that eta is a unitless parameter, right? So maybe viscosity/eta is not the right term for it?

Just of a minor note, particles.delta_gamma should save only the deviatoric part, that is gamma_dot_S. I realized that this is not explained anywhere... Please also consider renaming Delta_GAMMA to lowercase letters.

[atpellet]

Hi @larsblatny,

I provided a solution for the MCC_RMA_EXPLICIT() duplicate issue by passing delta_gamma as an optional pointer. It is not the most elegant or robust solution by this way no other change is needed.

They were some dimension issues with my viscosity, I have corrected and formatted it better. It now has a form often called Duvaut-Lion viscosity (there is extensive research on the subject). Basically, the viscoplastic strain is computed from the distance between the trail state and it's projection to the yield surface. The input parameter is the viscosity eta (Pa.s), but the compuations use a relaxation time defined as t_relax = eta/E (Youngs modulus).

I am not sure I fully understood your issue when saying One concern I have is in line 653 and 654 of plasticity.cpp where dydp and dydq (which you need to shift the (p,q)-point from the yield surface to a point away from the yield surafce given by eta) are calculated based on the point on the yield and not on the n+1 point.. Is it still a concern of yours ? If yes, can you develop ?

Cheers.
A. Pellet

PS : In my case, particles.delta_gamma stores the full plastic norm. I guess it is not an issue since this value is never used and just outputted in particles.ply ?

[larsblatny]

Hi @atpellet

Sorry for the late reply -- it took me a while to go through the math here:

The DL model is
$\sigma - P(\sigma) = t_{relax} \cdot (C : \dot{\varepsilon}^P) $
where $P(\sigma)$ is the projected state onto the yield surface. Typically $P(\sigma)$ is the closest point on the yield surface in principal stress space (and in this case the model would the equivalent to Perzyna with an associative flow rule). For this reason, let's denote $\sigma^{yield} = P(\sigma)$ in the following.

It is important to keep in mind that, in the equation above, $\sigma$ is not the trial stress state $\sigma^{trial}$. In fact, we have
$\sigma = \sigma^{trial} - C : \dot{\varepsilon}^P \cdot \Delta t$
Inserting this above and we obtain
$\sigma = a \cdot \sigma^{trial} + (1-a) \cdot \sigma^{yield}$ where we defined $a = \frac{1}{1+\Delta t / t_{relax}}$
or using the strain space and Matter notation,
$\varepsilon^{E, n+1} = a \cdot \varepsilon^{E, trial} + (1-a) \cdot \varepsilon^{E, yield}$

In comparison, the implemented model seems to be (correct me if I'm wrong):
$\varepsilon^{E, n+1} = \varepsilon^{E, trial} - \Delta \gamma \frac{\partial y}{\partial \sigma} \cdot \Delta t / t_{relax}$
I'm not sure how to interpret these variables so I will assume something similar (general) like this
$\varepsilon^{E, n+1} = \varepsilon^{E, trial} - b \cdot \dot{\gamma} \Delta t \frac{\partial y}{\partial \sigma} $
where $b$ is some dimensionless factor.
I'm also going to assume, here, which seems to be the case when looking at the code, that $\Delta \gamma = \dot{\gamma} \Delta t$ is the "strain step" between "trial" and "yield" (which does not define the plastic strain rate - this is defined between "trial" and "n+1"). If my math is correct, the implemented model would be a DL model if
$b = 1-a = \frac{1}{1+t_{relax}/\Delta t}$
assuming an associative flow rule / "closest point on yield in principal stress space".

If you could double-check my math/reasoning here and make the appropriate changes, that would be amazing! This would also include changing what you save as "delta_gamma".

Regarding the dydp and dydq, the evaluation of p and q seems to be done at "yield", which I suppose would make sense in the view of DL. There is the complication of the hardening law, but I think it is okay as it is.

Thanks a lot for all your work with this and let me know if anything is unclear. These are aspects of viscoplastic theory which I find super interesting!

[larsblatny]

With delta_gamma being the "strain difference between trial and n+1 found from the inviscid rma-function", and b = 1/(1+t_visc/dt), this is how I believe the code should be

T eps_dot_pl_vol_inst = - b * delta_gamma * dydp_yield;
T eps_dot_pl_dev_inst =   b * delta_gamma * dydq_yield * (1/d_prefac);

particles.eps_pl_vol[p] += eps_dot_pl_vol_inst;
particles.eps_pl_dev[p] += eps_dot_pl_dev_inst;

particles.delta_gamma[p] = d_prefac * eps_dot_pl_dev_inst / dt;

Note I used t_visc as this is already a parameter used for the viscous models and therefore we don't need to introduce another one :)

Please let me know if I made a mistake!

[larsblatny]

Hi @atpellet I reformulated the code according to my last comment, and I also removed the need for returning delta_gamma from the helper RMA function since what is needed can be calculated efficiently from p and q. Let me know if you have any comments, otherwise I'm fine with proceeding with merging.

[atpellet]

Hi Lars,

Sorry for not answering your comment.

Thank you for your inputs, it seems correct to me too. I did not answer earlier as there still are some details bugging me.

1. Strain decomposition :

The original formulation of Duvaut and Lion is
$$\dot\varepsilon = A\ :\dot\sigma + \frac{1}{2 \mu} [\sigma - P(\sigma)]$$
with $\dot\varepsilon$ the total strain.

Knowing that $A\ :\dot\sigma = \dot\varepsilon^E$, we rewrite this as
$$\dot\varepsilon^{vp} = \frac{1}{2 \mu} [\sigma - P(\sigma)]$$
by assuming decomposition of strains : $\varepsilon = \varepsilon^E + \varepsilon^{vp}$.

But this assumption is not true in finite strain theory. This bugged me aslwell in equation 2.67 of your thesis. When do/can we assume this decomposition in the framework of finte strain theory ?

2. Switching the formula

Assuming the previous relation, we have
$$\dot\varepsilon^{vp} = \frac{1}{2 \mu} [\sigma - P(\sigma)]$$

But the formula generally used (with no justification) written otherwise, as you did :
$$t_{relax} . (C\ :\dot\varepsilon^{vp}) = \sigma - P(\sigma)$$

By doing this, we assume that
$$2\mu . \dot\varepsilon^{vp} = t_{relax} . (C\ :\dot\varepsilon^{vp})$$
which instinctively seems correct. But we are comparing the scalar $2\mu$ and the scalar / 4th order tensor product $t_{relax} . C : $. Do you know of a formal way to switch from one expression to the other ?

Alexandre.

PS : Could you share if you based your reasoning on one particular reference ?

[atpellet]

Also,

I quickly reviewed your modifications. Everything seems fine. I just have some reservations regarding variable names.

• I think the relaxation time should be named as such. Here, you use the name perzyna_visc, which is very confusing. I think the variable should be renamed t_relax or at least t_visc in both models (Perzyna and Duvaut-Lion).

• The names of the models could also be changed for more clarity. MCCVisc and MCCVisc2 could become MCCPerzyna and MCCDL or MCCDuvautLion.

• Following our discussion, we could write tmp as b and write \\ tau = (1-b).tau_trial + b.tau_yield somewhere in the code in plasticity.cpp, for future readers.

If you agree with these suggestions, I can implement them ASAP.

Alexandre.

[larsblatny]

Hi @atpellet, excellent comments. I have now removed MCCVisc2 and merged it with MCCVisc, introducing the user option use_duvaut_lion_formulation. I renamed tmp, adding the comment you suggested. I agree that we should rename perzyna_exp and perzyna_visc. I suggest renaming them to viscous_exponent and viscous_time, but I will do this in another PR to avoid merging issues.

Regarding the "physics", I don't have a particular reference, I just derived it from the definition. There is for example a similar derivation here.

In finite strain theory, we need to "replace" the strain time derivative by the velocity gradient, and the elastic and plastic velocity gradients are additive. See Eq 2.53 in my thesis. Eq 2.67 is a one-dimensional special case.

Regarding the C : vs 2 * mu issue, I believe that the former is just a more general formulation of the latter. If we consider a VM or DP yield surface with a VM plastic potential (g=q), these two formulations will give the same (in the sense that we consider only the shear strain rate and shear stresses in the formulation). Do you agree?

[atpellet]

Hi @larsblatny, great solution to add the DL option in all models.

Regarding my derivation concerns, I agree it works in a shear-focused model, or any model reduced to a scalar relation.

I guess the idea of the DL model is just used in a convenient way in any framework without justification, thus we can use it aswell in the convenient formulation in our finite strain model. I guess it would be difficult to derive the same formula using the additivity of velocity gradients instead of strain time derivative ?

Anyway, I probably just have to assume the problem is well posed in this framework too, which feels awkward 🫤

[larsblatny]

Hi @atpellet
Yes, exactly. So the model would be something like $\sigma-P(\sigma) = t_{relax} C : l^P$. Or which formula are you thinking about deriving?
Well-posedness is tricky, even in the small strain case there are still unanswered questions here.