thelfer / MFrontGenericInterfaceSupport

This project aims at providing support for MFront generic behaviours. This project can be embedded in open-source and propriary sofware
37 stars 35 forks source link

Implementation of an interface to the general purpose multibody dynamics solver MBDyn #132

Closed octave-user closed 3 weeks ago

octave-user commented 5 months ago

Dear Mr. Thomas Helfer,

I would like to integrate your MFrontGenericInterfaceSupport library into the general purpose multibody dynamics solver MBDyn. So far, the small strain option is already working. However, the finite strain option seems to be more challenging. On one hand, I'm not 100% sure if I got the order of the tensor components right, and on the other hand I'm missing an option for the tangent operator. So, the derivative of the second Piola Kirchhoff Stress versus the deformation gradient "DS_DF" would be my favorite option.

You can find the current state of my implementation here: MFrontGenericInterfaceCSL and here: ComputeStressElasticDefGrad I'm also providing a link to my reference implementation of Mooney Rivlin which I'm using to validate the Signorini.mfront model. So, far I was able to run a simple uni-axial tensile test with this model, but the result was not as expected. Probably, I was using Cauchy stress instead of First Piola Kirchhoff stress, although the stress_measure was set to PK1. MooneyRivlinStressTensor

See also this link mbdyn-input-develop.pdf at page 313.

Any help would be appreciated.

Best regards, Reinhard

thelfer commented 5 months ago

@octave-user This is seems a very interesting project !

The order of the strain tensor and deformation gradient is given here: https://thelfer.github.io/tfel/web/tensors.html

The finite strain case shall not be that challenging. The main difference is in the load function which requires to specify the stress measure and the consistent tangent operator.

    auto o = FiniteStrainBehaviourOptions{};
    o.stress_measure = FiniteStrainBehaviourOptions::PK2;
    o.tangent_operator = FiniteStrainBehaviourOptions::DS_DEGL;
    const auto d = load(o, library, behaviour, h);

You can use the isStandardFiniteStrainBehaviour function to know if you shall use the "standard" load function or the one related to finite strain behaviour.

For example, this is the code that I wrote for the MFEM-MGIS project:

namespace mfem_mgis {

  std::unique_ptr<Behaviour> load(const std::string& l,
                                  const std::string& b,
                                  const Hypothesis h) {
    using namespace mgis::behaviour;
    if (isStandardFiniteStrainBehaviour(l, b)) {
      auto opts = FiniteStrainBehaviourOptions{};
      opts.stress_measure = FiniteStrainBehaviourOptions::PK1;
      opts.tangent_operator = FiniteStrainBehaviourOptions::DPK1_DF;
      return std::make_unique<Behaviour>(mgis::behaviour::load(opts, l, b, h));
    }
    return std::make_unique<Behaviour>(mgis::behaviour::load(l, b, h));
  }  // end of load

}  // end of namespace mfem_mgis

I won't have the time to have a look at your code this week, but if you want, we could have a short meeting next week.

octave-user commented 5 months ago

OK, the order of the tensor components and also the "load" function which I was using should be correct. Thank you!

But now I just saw, that the material parameters C01 and C10 were not consistent with my Mooney Rivlin parameters C1 and C2. So, the result seems to be correct now, but I need to check it in more detail.

thelfer commented 5 months ago

Great ! Let'us keep in touch and don't hesitate to make feedbacks or suggest proposals of improvement.

octave-user commented 5 months ago

Well, it would be nice to have an additional tangent operator, as I mentioned before. I think that such an option already exists, but it seems that it is not exposed by the generic interface. So, I would like to write something like this: o.stress_measure = FiniteStrainBehaviourOptions::PK2; o.tangent_operator = FiniteStrainBehaviourOptions::DPK2_DF; Actually it is not really necessary, but I could simplify my code and reduce some overhead.

octave-user commented 5 months ago

It seems that the First Piola Kirchhoff Stress Tensor PK1 is returned in a transposed format. Instead of T=(t11, t22, t33, t12, t21, t13, t31, t23, t32) I got T=(t11, t22, t33, t21, t12, t31, t13, t32, t23) Since this is not a symmetric tensor, the shear stress will be incorrect.

thelfer commented 5 months ago

Hi @octave-user, we use PK1 in FEniCS, XPER, Manta and many other projects. So far, nobody reported any issue. Could it be a confusion with the nominal stress ?

octave-user commented 5 months ago

Well, the difference may not be that obvious. I'm using a test case which was proposed by Konstantin Sedlan in 2000. I is about torsion tension coupling of an incompressible cylinder and a two parameter Mooney Rivlin material. That material is equivalent to MFront's Signorini model, provided that the parameter C20 is set to zero. I'm attaching a Maxima script which can be used, to generate the analytical expressions for the axial force and the torque. mooneyrivlin.tar.gz If you run that script, you get the following expressions for the torque M and the axial force N:


M=%pi*(C10+C01)*D*Ra^4

N=-(%pi*(C10+2*C01)*D^2*Ra^4)/3

If I'm using the Cauchy stress tensor, the output from MFront's Signorini model matches the analytical expression. However, if I'm using the First Piola Kirchhoff Stress Tensor and convert it back to Cauch stress, then the axial force is still correct, but the torque is no longer proportional to the twist D. But if I transpose the unsymmetric tensor, I get the correct result.

See the following test code: solidcsltest.cc

octave-user commented 5 months ago

See also page 102 of the following document: DissertationKonstantinSedlan.pdf

thelfer commented 5 months ago

If the Cauchy stress tensor is correct, then the only source of error would be convertCauchyStressToFirstPiolaKirchhoffStress from the TFEL/Math library which was generated with a CCAS system (maxima).

The code is a bit clumsy, but I checked a few terms and they seem correct.

It can also checked that those terms would be different if this function returned the transpose of the First Piola-Kirchhoff stress.

thelfer commented 5 months ago

What math library are you using ? Rare are math library whose first array index is 1...

thelfer commented 5 months ago

Are you sure of this line:

  return F * PK1 / F.dDet();

I have:

P = s F ^{-T} J => s = (P * F^{T})/ J

s being symmetric, we can transpose the last relation:

s = F * P^{T} / J

With this (hopefully correct) relation, you don't have to tranpose the MFront result.

octave-user commented 5 months ago

OK, I used the formulas from this book:

Richard Schwertassek Oskar Wallrapp Dynamik flexibler Mehrkörpersysteme Methoden der Mechanik zum rechnergestützten Entwurf und zur Analyse mechatronischer Systeme Friedr. Vieweg & Sohn Braunschweig/Wiesbaden image

thelfer commented 5 months ago

It seems that the definition of PK1 is mostly a matter of convention. Some authors used the definition that I gave earlier, but other used the transposed version. I don't know which notations are used by Schwertassek et al. but I guess that:

If so, their SIGMA is the transpose of "my" PK1.

octave-user commented 5 months ago

Probably you are right. Maybe it would be a good idea to add a comment in the source code or documentation (e.g. in the class FiniteStrainBehaviourOptions).

thelfer commented 5 months ago

To be honest, I would expect most calculations with MGIS in many solvers to fail if I made such a mistake :) But you are right about the documention. I'll do that ASAP. Thanks again for your feed-backs !

thelfer commented 5 months ago

I added a note on the convention used to define PK1. Let me know if you need more help

octave-user commented 5 months ago

OK, now it is clear. Thank you very much for your support.

My next step would be, to include a few MFront behaviors in my testsuite and to update the user manual including examples and references to your work. After that, I would like to ask you, to review the related section of the manual and to add also a reference to MBDyn at https://thelfer.github.io/mgis/web/index.html.

thelfer commented 5 months ago

Hi @octave-user, not problem. Tell me when it's ready. If you're satisfied with MFront and MGIS, don't hesitate to make a talk a the MFront user days !

octave-user commented 5 months ago

Please have a look at the following sections of the attached draft version of MBDyn's manual and let me know if you have any suggestions!

mbdyn-input-develop.pdf

thelfer commented 5 months ago

@octave-user You manual seems fine to me. I am a bit surprised that you have to specify a priori if the behaviour is finite or small strain, or the number of parameters.

octave-user commented 5 months ago

Well, MBDyn's constitutive laws and solid elements are based on templates with a fixed number of components of the stress/strain tensor (e.g. 6D for Green Lagrange Strain/Second Piola Kirchhoff Stress and 9D for Deformation Gradient and First Piola Kirchhoff Stress). Therefore, elements compiled for a 6D strain tensor cannot be used with MFront's finite strain behaviours which require a 9D strain tensor. The other way round would be possible somehow, but it is not implemented yet.

octave-user commented 5 months ago

In the meantime, I was playing a bit with examples from the MFrontGallery. Especially the SignoriniHyperviscoelasticity model seems to be very useful for my purposes. It worked without any issues, provided that all the viscoelastic parameters (g[0], g[1], g[2]) are set to zero. However, if I'm going to increase those parameters to a significant level, the solver is not able to converge any more. In contradiction to that, the GeneralizedMaxwell works without any issue, even if I raise the ViscoelasticShearModulus to very high values.

Does it mean that the viscoelastic contribution of the tangent operator for SignoriniHyperviscoelasticity is not consistent?

You can find a listing of my input file below: constitutive law: 1, name, "solid1", 9, mfront finite strain, library path, "/opt/mfront-gallery/lib/libHyperviscoelasticityBehaviours-generic.so", name, "SignoriniHyperviscoelasticity", properties, 10, "C10", 3.2307692307692310e+10, "C01", 8.0769230769230776e+09, "BulkModulus", 1.7499999999999997e+11, "C20", 0.0000000000000000e+00, "tau[0]", 1.0000000000000000e+02, "tau[1]", 1.0000000000000000e+03, "tau[2]", 1.0000000000000000e+04, "g[0]", 0.0000000000000000e+00, "g[1]", 0.0000000000000000e+00, "g[2]", 0.0000000000000000e+00;

thelfer commented 5 months ago

@octave-user Would you tell me what kind of loadings you are considering ?

octave-user commented 5 months ago

It is about torsion of a bar with rectangular cross section and one rbe3 element at each end. See the attached files (using the MooneyRivlin model): image solid.zip

thelfer commented 5 months ago

sorry, I don't have much time to learn the syntax of MBDyn's input files, but it seems that you make a time step of 1 (second ?). It this consistent with the internal times of the viscoelastic law ?

octave-user commented 5 months ago

OK, I think that I completely misunderstood the meaning of g. I thought, that it is something like the ViscoelasticShearModulus in the GeneralizedMaxwell model. But now I saw that it's unit should be rather non-dimensional. Since my values were scaled in [Pa], this was causing numerical truncation errors due to almost infinite stiffness. H[i] = e[i]*H[i]+g[i]*(1-e[i])/(dtr)*(Si-Si_1); However, if it is non-dimensional, it is not logical why it was declared as a stress. And also tau should be declared as time instead of stress, or am I wrong?.

@MaterialProperty stress tau[Nv];
@MaterialProperty stress g[Nv];
octave-user commented 5 months ago

Based on a finite difference approximation of the Jacobian, I can confirm that the tangent operator of the SignoriniHyperviscoelasticity model should be correct. I can also confirm, that the results from SignoriniHyperviscoelasticity are consistent with the GeneralizedMaxwell model for uni-axial tension and small strains.

So, the only remaining issue is about the misleading declaration of g and tau. See the following pull request: https://github.com/thelfer/MFrontGallery/pull/48

thelfer commented 5 months ago

Based on a finite difference approximation of the Jacobian, I can confirm that the tangent operator of the SignoriniHyperviscoelasticity model should be correct.

Great !

So, the only remaining issue is about the misleading declaration of g and tau.

I'll have a look at it. I will try to enable @UseQt true; to see if the compiler can detect correct unit usage.

octave-user commented 4 months ago

Dear Thomas Helfer,

I would like to inform you, that the new mfront interface was merged into the main branch of MBDyn. So, I think that you could add https://www.mbdyn.org/ to the list of projects using MGIS (https://thelfer.github.io/mgis/web/index.html).

Reinhard

thelfer commented 4 months ago

Great. Do you have:

so I could made some kind of announcement on LinkedIn (or if you do one, tell me, so I'll republish it)

octave-user commented 4 months ago

Well, you can use this one: MFront-Signorini-Hyperviscoelasticity.mpeg

octave-user commented 4 months ago

And here is a description from https://www.mbdyn.org/:

MBDyn is the first and one of the few full-featured free* general purpose Multibody Dynamics analysis software, released under GNU’s GPL 2.1 (get a cached copy here).

It has been developed at the Dipartimento di Scienze e Tecnologie Aerospaziali (formerly Dipartimento di Ingegneria Aerospaziale) of the University “Politecnico di Milano”, Italy.

MBDyn features the integrated multidisciplinary simulation of multibody, multiphysics systems, including nonlinear mechanics of rigid and flexible bodies (geometrically exact & composite-ready beam and shell finite elements, component mode synthesis elements, lumped elements) subjected to kinematic constraints, along with smart materials, electric networks, active control, hydraulic networks, and essential fixed-wing and rotorcraft aerodynamics.

thelfer commented 4 months ago

Many thanks. See this post: https://www.linkedin.com/feed/update/urn:li:share:7198382956053176321/

thelfer commented 2 months ago

Hi @octave-user, Could we close this issue ?