TFEL/MFront + FEniCS - Cosserat elastoplastic behaviour

[tamaradanceva]

Hi,

I am a Ph.D. student at the Basque Center for Applied Mathematics in Bilbao, Spain. I am also part of the ENABLE project, an Innovative Training Network (https://enable-project.com). My specialization is in High-Performance Computing, in the context of setting up an HPC environment and developing efficient methods for advanced multi-scale material behavior law models using strain-gradient plasticity models, in the open-source Finite Element platform, FEniCS.

I have been exploring the possibility of using MFront to implement a Cosserat elastoplasticity model in the framework of small deformations as a starting point. The next objective, ideally, would be implementing a similar Cosserat model for large deformation. Now, I have seen that I can run simulations with multiple threads. My question is if there is a workflow (experimental or otherwise) that would enable me to run in parallel? I am following the progress on Dolfin-X. However, in the meantime, is there any alternative way you are aware of to accomplish my goal? Implementing a solid mechanics library like FEniCS Solid Mechanics, from scratch is not feasible since this is only one of a series of tasks that I need to do within my project (with a duration of 3 years).

Your work has been of great help to me so far. I would be very grateful for any input that you can give me.

Thanks,
Tamara

[thelfer]

Dear Tamara,

Thanks for your interest. I think that we can do something interesting around the mgis.fenics python module which has mostly been developed by @bleyerj. See this page for details: https://thelfer.github.io/mgis/web/mgis_fenics.html

For Cosserat elastoplasticty, I need to have a closer look at the behaviour that you want to implement to tell you if you can implement it right now or if some features are missing in TFEL/MFront. Have you any reference illustrating what you want to achieve ? Anyway, as I am interested by such developments, I will implement them if needed.

Performance wise, we currently rely on the MPI parallelisation provided by FEniCS and on each MPI process, we can use the multithreading parallelisation provided by MGIS. But there is certainly room for improvements. This is something that we will tackle when porting the module on dolfin-X. We already have some ongoing discussions with a member of the FEniCS steering counsil working on dolfin-X. I also plan to start working of porting MFront and MGIS to GPUs but this is probably a matter of two to three years before having something solid on this topic.

Concerning the FEniCS Solid Mechanics project, there is an experimental attempt in MGIS to write something similar (some of the code is directly copied/inspired for FEniCS Solid Mechanics), i.e. interface MGIS and FEniCS at the C++ level. I must admit that I currently consider it as a dead end. Probably more efficient than mgis.fencis (we did not make any benchmark so it is a purely a priori statement), but almost all what makes the power of mgis.fenics module, i.e. its astonishing flexibility and ability to define very complex mechanical models, is lost.

Do not hesitate if you have further questions, we will be happy to help.

Regards,

Thomas

[bleyerj]

Hi @tamaradanceva ,
If I am not mistaken Cosserat models (elastoplasticity in particular) should be straightforward to implement using TFEL/MFront and the mgis.fenics module since generalized stresses are just a non-symmetric 2nd-order stress tensor and a non-symmetric 2nd-order couple stress tensor if I am not mistaken. This should be already supported.

If you give us the exact small strain model you want to solve we can try an implementation. I am not familiar with large strain Cosserat models but we can also have a look if small strain is OK.

Strain-gradient models should be harder for two reasons (again if I am not mistaken):

• generalized strains involve the strain gradient which is a 3rd-order tensor, these kind of objects are not currently supported and should be represented in vectorial notation but there might be some missing features in TFEL/MFront and mgis.fenics. Again, an example of what you would like to do would help us in identifying these missing features
• regarding finite-element interpolations, C1 continuity is required. Such elements are however not currently supported and you would need to resort to some relaxation or a Discontinuous Galerkin approach. Not impossible but more work to do.

[tamaradanceva]

Hi, @thelfer and @bleyerj,

Thank you very much for the answers. I will get back to you if I have any more questions since I will try implementing the models. First I am trying to run some tests in parallel.

Here is a link to a preprint of an article by my colleague in the modeling work package that is going to be published in the next issue of the Journal Continuum Mechanics and Thermodynamics: https://zenodo.org/record/4059439. It describes in detail the derivation of an Elasto-Visco-Plastic Cosserat model in small deformation applied on the Hat-Shaped Specimen benchmark (a loading that causes the localization of the shear deformation in a band - called the Adiabatic Shear Band). In the Cosserat model, the rotation of the microstructure is a vector, and the Cosserat wryness tensor (equivalent to the micro-strain gradient tensor) results in a non-symmetric 2nd-order tensor. My colleague decided to use the Cosserat instead of the strain-gradient approach (that figures in the project proposal) since it is more robust and less approximative, in short, in his words. It can mimic the same properties and finally, the strain gradient can be obtained with a choice of parameters.

Best,
Tamara

[RaffaRus]

Hi @thelfer and @bleyerj ,

First of all many thanks for the support!

I am collaborating with @tamaradanceva on the implementation of the elasto-plastic version of the Cosserat media under Small Deformation. We are in the same European Project ENABLE, and I am currently implementing the Large Deformation version (Hyper-Elastic, Visco and Elasto Plastic) in Z-Set (ZEBULON).

I am not familiar with the syntax by which the *.mfront files are written, but I recognize that many of the *.cxx structure follows similar logic than in Z-set. The main issue that we are facing is the selection of an appropriate class of material behavior that would allow us to define additional Grad, additional Flux and additional material tangent matrix. Here I report a short summary of the kinematic and material laws governing our problem:

// Kinematic
u = displacement field defined in R3;
w = micro-rotation field defined in R3;
The strain is e = grad(u) + permutation(w);
The wryness(curvature+torsion) is k = grad(w);

// Elasticity
sigma and m are the stress and "couple stress";
sigma = lambda*tr(e)*I+2*mu*sym(e)+2*mu_c*skew(e);
m = alpha*tr(k)*I+2*beta*sym(k)+2*gamma*skew(k);
E and C are the fourth order tensors, elastic operators, respectively related to strain/stress and wryness/couple stress;

//Plasticity
s_dev = deviatoric stress;
p is the palstic multiplier;
sigma_eq = sqrt(a_1*s_dev:s_dev + a_2*s_dev:transpose(s_dev) +  b_1*m:m + b_2*m:transpose(m));
f  = sigma_eq - R;
R = H*p;
                       N*E*dot(e) + N_c*C*dot(k)
     dot(p) = ---------------------------------  ;
                        H + N*E*N + N_c*C*N_c
dot(e)^p = dot(p) * (a_1*s_dev + a_2*transpose(s_dev))/(sigma_eq);
dot(k)^p = dot(p) * (b_1*m + b_2*transpose(m))/(sigma_eq);

We unsuccessfully tried to use the DefaultDSL and DefaultGenericBehaviourDSL material classes but got into issues when trying to define additional material tangent matrices.

Many thanks in advance for the help.

Best,
Raffaele.

[thelfer]

Hi @RaffaRus,
I am currently out of office. I’ll have a deeper look at it next week. Some tensorial operations seem missing right now, but nothing critical. The DefaultGenericBehaviourDSL shall do the word, but it is fairly undocumented right now. Would you be ok to help us write a new page for the MFront gallery or a new tutorial for mgis.fenics ?

[RaffaRus]

Thanks a lot for the assistance! Yes, if I can be of any help writing some tutorial I would be glad to do so. Let us talk next week then.

Have a nice weekend.

Best,
Raffaele.

[bleyerj]

Hi @RaffaRus, I think that this could be feasible right now. Could you share with us your best attempt and we can have a look and get back to you.
Best

[RaffaRus]

Hi @bleyerj , thanks for your help! In the attachemnts our attempt to model the behavior using the DefaultGenericBehaviour class. There are few comments related to the issues we are facing. Please let us know if they are not clear.

Best.

CosseratIsotropicLinearHardeningPlasticity.mfront.txt

[bleyerj]

Indeed, there seems to be only partial support of Tensor/Tensor tangent blocks in MFront which makes straightforward implementation of Cosserat EP behaviour not possible right now. This should be fixed soon (I opened the following ticket https://sourceforge.net/p/tfel/tickets/252/)

[thelfer]

Hi
There are indeed a few features missing. The tensorial objects exists, but not associated views in externally allocated memory. I’ll Check this out and let you inform.

[thelfer]

@RaffaRus, @bleyerj. With a one line fix, here is a minimal example.

@DSL DefaultGenericBehaviour;
@Behaviour CosseratIsotropicLinearHardeningPlasticity2;
@Author Raffaele Russo, Tamara Dancheva;
@Date 26 / 10 / 2020;

@Gradient DeformationGradientTensor F;
@Flux StressTensor PK1;
@Gradient DeformationGradientTensor kappa;
@Flux StressTensor m;

@Integrator {
dPK1_ddF = t2tot2<N,real>::Id();
dm_ddkappa = t2tot2<N,real>::Id();
}

The fix is available shortly on github.

[thelfer]

It would be nice to detail the algorithm. I understand that this is a fully Euler backward algorithm, is it ?
If this is the case, then the DefaultGenericBehaviour DSL is already sufficient (no work on the MFront side). However, some functions shall be added to the TFEL/Math library (for example, the skew function and its derivative).

@RaffaRus would you start writing a small document describing the principle of virtual box and the behaviour integration. For the better would be to write it in pandoc (see for example the docs/web/ExplicitBerveillerZaouiPolyCrystals.md file in the TFEL sources), but LaTeX is quite fine, so that we can discuss the better choice for the definition of the gradients and thermodynamic forces, tangent operator blocks.

[bleyerj]

hi, here is a first proposal which compiles with the current master. I didn't check it yet though

CosseratIsotropicLinearHardeningPlasticity.mfront.txt

[thelfer]

Hi @bleyerj,

This is becoming nice but:

1. @PredictionOperator is not called before @Integrator, so ∂σ∕∂Δεᵗᵒ and ∂μ∕∂Δκ are not defined when you use them
2. there's no tangent operator. For this, I just want to know if there are operators missing.

[bleyerj]

yes, I was just sending this quickly, I had my kids around... I'll get back to it this week

[thelfer]

@bleyerj, come mine too. Also, was not sure if the proposed algorithm is fully implicit or fully explicit... That's why I proposed to start write a small paper to fix the equations (and the computation of the consistent tangent operator)

[tamaradanceva]

Hi @thelfer @bleyerj,

Thank you very much for your help so far. We have tested the version sent by @bleyerj and it doesn't compile, with or without the definition of the tangent operator. It results in the same error as before, from the check of the types of the gradient and the flux when creating the tangent operator. This is including with the latest version.

BehaviourDSLCommon::getBehaviourConstructorsInitializers:tangent operator blocks associated with the derivative of 'σ' (of type 'Tensor') with respect to 'εᵗᵒ' (of type 'Tensor') is not supported

This is just to check that we are on the same page. Meanwhile, @RaffaRus is going to provide a small document with the description of the full model, the virtual work and the behaviour integration.

[thelfer]

Hi @tamaradanceva, @RaffaRus
Indeed, @bleyerj did a great job. I think that we are close to something working.
Have you updated your TFEL version. @bleyerj and I used the latest commit of the master branch.

[tamaradanceva]

Hi @thelfer and @bleyerj ,
Sorry, my mistake, I am on the latest master now and it compiles! There are some more changes that we need to do on the FEniCS side and a part of the plastic tangent operator. We will test it further and let you know of the progress.

[bleyerj]

OK, let me know if you need some help on the FEniCS side, in particular it may be possible that the Cosserat model needs some reduced integration strategy to prevent any numerical locking...

[bleyerj]

@bleyerj, come mine too. Also, was not sure if the proposed algorithm is fully implicit or fully explicit... That's why I proposed to start write a small paper to fix the equations (and the computation of the consistent tangent operator)

@thelfer the algorithm will be implicit but with a closed-form analytical expression as for the standard von Mises with isotropic linear hardening. The are however alternative models for Cosserat plasticity which are more like a multisurface plasticity e.g. a von Mises surface in the stress-space intersected by a von Mises-like surface in couple-stress space. Of course, you can also have more complex yield surface shapes like Drucker-Prager, etc...

[thelfer]

@thelfer the algorithm will be implicit but with a closed-form analytical expression as for the standard von Mises with isotropic linear hardening.

Ok, I just did not take the time to write it properly. So your implementation already shows that the current DefaultGenericBehaviour DSL is enough, one just have to work on the tangent operator side. Good.

For the moment, the ImplicitGenericBehaviour DSL don't accept (non symmetric) tensor variables as integration variables, so I feel a bit relieved that the DefaultGenericBehaviour DSL does the job. Extending the ImplicitGenericBehaviour DSL would have been cumbersome. This is because the increments are views on part of a great array of values which is manipulated by the Newton algorithm and the derivatives are views of a big jacobian matrix. Such views don't exists right now for tensor, t2tost2, st2tot2 and t2tot2 objects. I plan to implement them in TFEL-4.0 (the next version) which will be rewritten in C++-17. Backward compatibility will be guaranteed: every MFront file will compile in this new version, but internally everything will be overhauled for greater efficiency and maintainability (and probably support of GPUs).

[RaffaRus]

Hi @thelfer and @bleyerj ,

Here I attach a small report on the model we are trying to implement. We are also close to something working thanks to your help. We are going to let you know asap. Please let me know if I missed some important information in the report.

Cheers,

Raffa.
Cosserat SD - Model.pdf

[tamaradanceva]

Hi @thelfer, @bleyerj

Here is a link to the case we are working on, https://github.com/RaffaRus/Cosserat-Fenics/tree/master/Tensile_Ti-6Al-4V-Cosserat. We did some benchmarking on a tensile test with classical elastoplasticity and we are trying to extend the same example with Cosserat for a start. It should not make a difference if it is done correctly. I want to ask you specifically about the transfer of variables between MFront and FEniCS (gradients, thermodynamic_forces, and K). I think that as it is, it is not correct and the displacement rotation solution is diverging partly because of that and other potential issues. From our understanding and looking into the code, for example in the case of K, which I suppose it is the same for the rest, the values alternate Ct[0][0], Dt[0][0], C[0][1], D[0][1].. Is this is a correct supposition?

Best,
Tamara

[bleyerj]

dear @tamaradanceva and @RaffaRus , as discussed today with @thelfer, it seems that :

• you did not consider coupling tangent terms dsig/dkappa and dmu/deps which should be considered for the consistent tangent operator (note that such coupling is activated only in the plastic case)
• contrary to what is obtained for classical J2 plasticity with isotropic elastic behaviour, the implicit system cannot be solved analytically in the Cosserat case. Indeed, in classical J2 plasticity one has n_{elas}=n i.e. the unknown normal n is equal to that obtained with the elastic trial state. This does not seem to be the case with the Cosserat behaviour. Hence, as feared initially, we should resort to the ImplicitGenericBehaviour DSL which will not be available with non-symmetric tensors before TFEL-4.0
• finally, we suggest using the mgis.fenics project to avoid the troubles you just mentioned above. Here is an example of such a script that I wrote (note that it does not converge when plasticity is activated due to the issue mentioned above)

cosserat_mgis.zip

[thelfer]

Hi @tamaradanceva, @RaffaRus, @bleyerj

contrary to what is obtained for classical J2 plasticity with isotropic elastic behaviour, the implicit system cannot be solved analytically in the Cosserat case. Indeed, in classical J2 plasticity one has n_{elas}=n i.e. the unknown normal n is equal to that obtained with the elastic trial state. This does not seem to be the case with the Cosserat behaviour. Hence, as feared initially, we should resort to the ImplicitGenericBehaviour DSL which will not be available with non-symmetric tensors before TFEL-4.0

Thanks for checking this, I indeed came to the same conclusion, i.e. that no closed form solution exists for the implicit scheme in the Cosserat case. So we'll have to rely on the ImplicitGenericBehaviour DSL. I will do the required developments after the release of TFEL-3.4 (end of month) in the master branch (i.e. the future TFEL-4.0 version) at the beginning of december.

In the meantime, we would already use an explicit integration scheme (based on the consistency condition deveopped in the document provided by @RaffaRus), but the consistent tangent operator will not be available. As a result, the global convergence will be only linear and not quadratic. @bleyerj Does FEniCS provide some kind of quasi-newton algorithm (BFGS ?) If not, we could try the Anderson acceleration scheme provided by TFEL which is very efficient (a variant is used in Castem for decades and we used this algorithm in various FFT solvers). What do you think ?

[thelfer]

@RaffaRus, @tamaradanceva, @bleyerj Here is an explicit implementation based on Equations 1.59, 1.60 and 1.61. More precisely, I adapted Jérémy's implementation.

Only the elastic operator is returned. We could implement the tangent operator given by Equations 1.57 and 1.58 (completed by the coupling blocks), however this is probably a bad idea: in standard elasticity, it is often advocated in the litterature that such tangent operator leads to a spurious oscillations of the residual. We shall better stick with the elastic operator and acceleration algorithms.

I did not test it so far, I'll see when I'll have some time. @bleyerj do not hesitate to give a try in the meantime.

Note that you will need the very latest version of TFEL to compile this behaviour, as I had to fix Issue #255 (https://sourceforge.net/p/tfel/tickets/255/).

ExplicitCosseratIsotropicLinearHardeningPlasticity.zip

P.S.: @bleyerj In your implementation, beware that (3/2) is equal to 1 in C++.

[bleyerj]

@tamaradanceva @RaffaRus
Thanks to @thelfer explicit implementation, I have now a working example using mgis.fenics with the SNES solver using Newton method (with the elastic operator) and a line search procedure
If you could check that with an analytical reference solution that would be great.
You will also need MGIS latest version, as I had to fix something to make it work

explicit_cosserat_mgis_fenics.zip

[thelfer]

@bleyerj Great work !

@RaffaRus, @tamaradanceva don't forget to update both TFEL and MGIS to make it work. At this stage, performances can still be improved (next month), but early tests can be made.

[tamaradanceva]

Hi @bleyerj, @thelfer

Thank you very much for your help! We will test this approach with giide and bending for which Raffa has a semi-analytical solution and let you know.

[thelfer]

@RaffaRus, @tamaradanceva, @bleyerj I just wanted to point this issue: https://sourceforge.net/p/tfel/tickets/256/
It may affect your computation. I solved it in the current version, I would recommend that you update your version.

[tamaradanceva]

@thelfer Thanks for the heads up, I updated it.

[thelfer]

The github has been updated at 20h43.

[tamaradanceva]

Hi @thelfer, @bleyerj

Here is a link to the glide test, https://github.com/RaffaRus/Cosserat-Fenics/tree/master/Glide-Cosserat-ElastoPlasticity with results from FEniCS and Zset, and a comparison with the analytical solution.

We are going to continue with the bending test, performance analysis, extending the model etc.

Thanks a lot for your help again!

[thelfer]

@tamaradanceva Thanks for the feed-back, that's quite nice.

As discussed before, performances are yet to be optimised and we expect two main improvements:

• by the end of the year (TFEL-4.0), an implicit implementation of your behaviour with a consistent tangent operator
• during next year, port of mgis.fenics to dolfin-x

Could you also have a look at the number of iterations to find the equilibrium in each cases ?

[tamaradanceva]

@thelfer,

About the improvements, I have them in mind. I was about to ask you about the plans for moving to dolfin-x, it is good to know. We just have an internal deadline at the beginning of December, for which, I would have to do an initial evaluation.

In the FEniCS+MFront case, we have [1,1,1,1,1,1,3,7,8,8,9,10,10,11,12,13,15,16,18] iterations. I will get back to you on Zset, I have to check with Raffa.

[thelfer]

@tamaradanceva That's good to have those limitations in mind (and the fact that they will hopefully be overcome in mid-term).
It seems that you don't use mgis.fenics yet. Your input files will be much shorter and easier to read/maintain/adapt, but that would also help use other FEniCS built-in solvers (quasi-newton, BFGS, etc...). @bleyerj don't hesistate to comment, you're far more experienced than I am on that topic.
@tamaradanceva @RaffaRus @bleyerj After the MFront User Meeting (Novermber the 25th), it would be nice if we could organize a small meeting about your project and discussing current results and future plans, don't you think ?

[tamaradanceva]

@thelfer For the glide we do, but for the rest not yet. We started exploring the options initially with the elasto plasticity demo, a tensile test, later realized it can be simplified so much with mgis.fenics .

Could we also attend the MFront User Meeting? I looked it up but found no mention. About organizing a meeting afterwards, absolutely.

[thelfer]

@tamaradanceva, @RaffaRus Of course you can attend the meeting. Just register at [email protected].
The agenda is available here: https://github.com/thelfer/tfel-doc/blob/master/MFrontUserDays/SixthUserDays/talks.pdf
The annoucement has been made on a dedicated mailing list of the project (I can add you if you want), the Twitter account of TFEL/MFront (https://twitter.com/TFEL_MFront) and my personnal LinkedIn account.

[tamaradanceva]

Hi @thelfer, @bleyerj

Firstly, firstly, I said I would get back to you on the number of iterations in Zset for the glide test. It is: [1,1,1,1,1,1,4,5,6,6,6,7,7,8,9,10,10,9].

Secondly, about the meeting, we propose next week Thursday (03.12), after 13:00, and Friday (04.12), any time you are available. Let us know if any of these slots would work for you or other slots that work for you.