Neo-Hooke implementation

[bengt]

Hi, dear Elmer friends,
since there have already been several posts concerning non-linear elastic material laws (e.g. for rubber) and I'd pretty much like to use such a model, too, I've gathered some information about this topic.
I'd like to share some ideas, hoping to launch an interesting discussion...
First of all, I think the current "non-linear" (or "large displacement") elasticity implementation within "ElasticSolve" refers to the "generalized" Hookean law being derived from the elastic strain-energy function, e.g. as reproduced in the attached pdf from FEBio theory manual (University of Utah, http://mrl.sci.utah.edu/febio-overview). Now, the second Piola-Kirchhoff stress can be derived from this function and, as far as I understand, this formulation is used in Elmer, too, if we take a look at the code:

Stress2 = 2.0D0*Lame2*Strain + Lame1*TRACE(Strain,dim)*Identity

(a commented version of the source is attached, too).
Using Neo-Hooke as a first non-linear material approach (see section 5.3 of pdf), would it be possible to slightly modify the current ElasticSolve in order to get a dedicated - let's say - ElasticSolverNeoHooke? As far as I understand, the Lamé parameters from linear elasticity could still be used, while only (??) the stress formulation would have to be modified . A very nice illustration of what we might need is given by François Labelle on http://wismuth.com/elas/elasticity.html:

S = lambda*J*(J-1)*C^-1 + mu*(I - C^-1), where

J = det(deformation gradient)
C = F^T * F = right Cauchy-Green deformation tensor
lambda, mu = Lamé constants

I've marked the according lines and symbols in the attached source file.
Trying to modify the stress formulation, however I just wasn't able to find out how "invertMatrix" is called (and used) correctly...
And...first of all: do you think this kind of implementation is correct and could actually be done without any major additional coding efforts??
I'd greatly appreciate your suggestions!
Kind regards, enjoy the weekend
Bengt

[hamed]

Hi Bengt
The material model implemented in Elmer is the Saint Venant-Kirchhoff model (SEE e.g. chapter 5 of "Finite Elements for Nonlinear Continua & Structures, Belytschko" ). However, I should confess that there are many lines in the source code of subroutine "StressCompose" that I can't understand. So I think it would be a very nice idea to discuss this solver in details.
Regards, Hamed

[bengt]

Hi Hamed,
thanks for your feedback! From what I read in the meantime I think it could be possible to modify the solver as discussed, but still don't know how to code the "invertMatrix" correctly...
Best regards, Bengt

[bengt]

Dear Peter and Juha,

do you think the Neo-Hooke implementation could be done as discussed above? I mean, wouldn't it be sufficient to modify the "stress2 = ..." - lines according to

S = lambda*J*(J-1)*C^-1 + mu*(I - C^-1)

in order to get a dedicated "NeoHooke-Solver"?
And if yes, would it be possible for you to give us some guidance on how to code the modified stress equation, especially on the correct use of "invertMatrix"?
I'd greatly appreciate your suggestions.
Best regards, Bengt

[mika]

Hi,

By exploring the code and taking the risk that I am not right, I would say that the finite elasticity solver treats the equations in the form where the constitutive equation describing the Piola-Kirchhoff stress is written as S = F G, with F the deformation gradient and G a response function defined in terms of the Lagrangian strain tensor E as

G = lambda*trace(E) I + 2*mu*E.

Having the constitutive equation in the form S = F G(C) generally relates to the requirement that the response is independent of the observer. It seems to me that a more standard and still observer-independent way to write the constitutive equation for describing the behavior of a hyperelastic material would, however, be to express the constitutive equation simply as

S = lambda*trace(E) I + 2*mu*E.

Considering the implementation of other constitutive models, I would start from differentiating the constitutive relation S = S(F) in order to find a feasible linearization strategy and then explore in more detail what modifications to the existing code are required. That is, I expect that the modifications needed depend essentially on the choice of linearization.

Best regards,
Mika

[bengt]

Hi Mika,

thanks so much for having a look at the code!! To be honest, I'm definitely not an expert in computational solid mechanics, so discussing things is rather difficult for me. Regarding the fem implementation of hyperelastic material behaviour I recently discovered the following online ressource:
http://solidmechanics.org/text/Chapter8 ... ter8_4.htm
As said, I'm not able to discuss things in detail, but from your own personal experience, do you think that the strategy being described there (for neo-hooke) could be used with the elmer code?
Trying to understand the ElasticSolver code, I've been partly able to understand the use of the different subroutines. However, would it be possible for you to briefly explain what "BackStressGeneralCompose" is doing in addition to "StressCompose" ("StressGeneralCompose" is oviously needed only in the case of non-Cartesian coordinate system)?
Again, thank you very much for your support!
Best regards, Bengt

[mika]

Hi,

First of all I would like to notify you that the documentation of the finite elasticity solver has eventually been done. I am attaching the relevant pages of the revised model documentation that concern finite elasticity. I hope that this limited piece of information helps in understanding what are the equations and how they are treated. The source code has already been revised today.

In addition to handling the acceleration term and body forces, the routine StressCompose treats the terms that arise when the Newton approximation (6.5) is substituted into the weak formulation. As you can see, the first three terms in the last line of (6.3) depend only on the previous iterate and thus contribute into the right-hand side vector. The remaining terms have an effect on the stiffness matrix computation. They are however implemented in the form that is based on using the next equation below (6.3). The routine BackStressGeneralCompose is never called, so I have not looked into it.

As anticipated, the constitutive equation has indeed the form where the first Piola-Kirchhoff stress is expresses as S = F G, with G giving the second Piola-Kirchhoff stress. If you now consider implementing some other constitutive model with a different definition of G, it seems to me that you could utilize the structure of the existing code provided you can evaluate the derivative DG(u)[v] in a closed form, as the equation (6.3) basically gives a linearization for any response function G.

Best regards,
Mika

[mika]

Hi,

The equation number (6.3) in my previous message should of course be (6.5).

-Mika

[bengt]

Hi Mika,

sorry for the late reply, I wasn't able to follow the forum during the last days. Again thank you so much for all your efforts regarding the documentation! Unfortunately I feel that I won't be able to do the rquired coding myself, so I'd like to ask all the interested forum members if someone already knows how the implementation could be done and if he or she would like to share ideas or results.
Kind regards and you all have a good week!
Bengt

[bengt]

Hi, dear Elmer team,

well, indeed I haven't been able to solve the neo-hooke material implementation my own and unfortunately we didn't get feedback from the community either.
Do you @Elmer think it might be possible to help with the implementation, or would that include complex modifications of the existing code? If it means major coding efforts, would it be then possible to include this feature in the list of official feature requests? I guess it would be of great help for many Elmer users to have such a material model available.
Best regards
Bengt