Holzapfel-Gasser-Ogden (HGO) material model - unable to reproduce Abaqus results

[zigadonik]

Hi,

As the title states, I am trying to reproduce a sample case for HGO in FEBio that is available in Abaqus documentation (INP file). I modified the input file (attached as feb) to load the whole upper surface with a displacement load (instead of a concentrated force and linear constraint as in the original input file) to make FEBio implementation easier (I did not find linear constraints in FEBio Studio, though I know they can be added to FEB file). I used the uncoupled HGO (incompressible) in FEBio, with the material parameters from Abaqus documentation (except for the shear modulus: c[feb] = 2 * C10[abq]) and changed the material axes to be able to use the same angle (49.98°). I can not get similar results (stress, reaction force in N1). What am I missing?

Versions: FEBio Studio 1.6.1 (FEBio 3.5.1), Abaqus 2021 SE.

I appreciate any help you can provide.

Žiga

[ateshian]

Hi Žiga,

Your analysis ran correctly. To get the same results as Abaqus, you should load the plot file (.xplt), select all the nodes on the right x-face of the bar, then select the X-reaction force as the plot variable to display, then go to Post->Integrate to get the total X-reaction force on this face. A window will pop up, showing the net reaction force as a function of time, let's call it F. Copy this curve to the clipboard and paste it into a spreadsheet.

You can also select one of the nodes on that same face and display its X-displacement using Post->New Graph...->displacement->X-displacement. Copy and paste this displacement onto the spreadsheet, let's call it d.

Since your model is symmetric about the X-plane, you need to create a new column in your spreadsheet which is equal to 2*d. Since the force F acts on only a quarter of the model (since the model was split in half along the Y- and Z-planes) you need to multiply this force by 4 (actually -4 to get a positive force), thus 4*F. If you now plot 4*F (ordinate) versus 2*d (abscissa) you get the same result as Abaqus (and as the Gasser-Ogden-Holzapfer 2006 paper from which the Abaqus graphs is copied).

Best,

Gerard

[zigadonik]

Hi Gerard,

Thank you for the swift reply and the clarification. I managed to get the same force-displacement curve as in the 2006 paper. I was a little puzzled by the difference in maximum stress in X-direction, where there was a significant difference. Abaqus requires C3D8H (hybrid) elements for the anisotropic hyperelastic model, whereas FEBio uses a different element formulation. Could element formulation be the reason for the difference in stresses?

Best,

Žiga

[ateshian]

Hi Žiga,

Yes, using different element types can lead to differences in the stress. For example, you can switch your model from hex8 to hex20 elements and you will see that the stress distribution changes, especially for the values of the minumum and maximum values of the normal stress along the X-direction.

Best,

Gerard

[athpa]

Hi all,

It is really helpful this post for us who working with HGO model. However I have some concerns about the possibility to reproduce the case which is available in Abaqus and especially the axial sample. As you can see in the attached file, after the sum of the reaction forces (x4 as ateshian mentioned) is 0.121941 which is far way of the curve that published in the paper. And also for the circumferential sample the curve seems the same but the forces not.

Did I do something wrong or miss something ?

In parallel I tried to apply force or rigid force to the model without success getting every time error of Negative Jacobian.

I would appreciate any help or remark on this topic.

Best,

Athpa

[ateshian]

Hi Athpa,

I tried running the first file (hgodisp1axial.feb) but there were a number of errors in the way the boundary conditions were prescribed. I tried fixing them and got the model to run (hgodisp1axial.fs2) , but I don't know that I did it right. Would you please share the exact boundary conditions and results from the ABAQUS model you are referring to? This will help me answer your questions.

Best,

Gerard

[athpa]

Hi Athpa,

I tried running the first file (hgodisp1axial.feb) but there were a number of errors in the way the boundary conditions were prescribed. I tried fixing them and got the model to run ([ATTACH]n22660[/ATTACH]) , but I don't know that I did it right. Would you please share the exact boundary conditions and results from the ABAQUS model you are referring to? This will help me answer your questions.

Best,

Gerard

Hi Gerard,

Thank you for your reply .

I upload the files in .feb because I can't upload txt or zip files.

Here is the linke https://abaqus-docs.mit.edu/2017/Eng...yper-load-disp where you can find the results of the published papers of Gasser et al (2006) which I am referring to and below the load/displacement curve that I am interested to.

image.png



Best,
Athpa

[ateshian]

Hi Athpa,

Thanks for sharing that web URL. I downloaded the ABAQUS *.inp files and imported the mesh from one of them into FEBioStudio and set up the analysis for the four cases ( adventitia_axial_kp226.fs2 adventitia_axial_k0.fs2 adventitia_circ_kp226.fs2 adventitia_circ_k0.fs2 ). When I plotted the FEBio results against those of ABAQUS I got the following agreement:
FEBio-ABAQUS-HGO-model.png
What you need to keep in mind is the following: The ABAQUS mesh uses symmetry about the x-, y- and z-planes. So, when you convert the results from the finite element analysis to the above plot you need to multiply the x-displacement of the moving face by 2, but the force by 4.

Another thing to keep in mind is that the case k=0.266 in the figure is really k=0.226 (as can be determined from the ABAQUS input files). The label 0.266 was used here for consistency with the typographical error appearing in the Gasser, Ogden, Holzapfel (2006) paper.

[ateshian]

The agreement between FEBio and ABAQUS for the case k=0.226 is excellent. The relatively small discrepancy in the response for k=0 is, in my opinion, a problem with ABAQUS. The reason I say this is because the case k=0 is equivalent to having two single fiber families. In FEBio, this can be reproduced using an uncoupled solid mixture with a Mooney-Rivlin ground matrix and two fiber-exp-pow-uncoupled fibers ( adventitia_axial_k0_mixture.fs2 adventitia_circ_k0_mixture.fs2 ). The results for the mixture analysis are identical to the results for the HGO model with k=0, so that gives me confidence that our implementation of HGO is correct.

To enforce strict incompressibility, I used the three-field method in the attached FEBioStudio files. So, if you plot the relative volume J as a function of deformation, you will find that J remains equal to 1 (within at least three decimal places) in the FEBio analyses. I don't know if ABAQUS enforces incompressibility that well, so it wouldn't hurt to check that on your end if you can.

As a final note, the ABAQUS solution does not agree with the corresponding figure published in the Gasser, Ogden, Holzapfel (2006) paper any better than the agreement between ABAQUS and FEBio shown above. Here is a plot of the ABAQUS results superposed over the GOH paper:
ABAQUS-GOH-paper.png

The discrepancy is greatest for the axial speciment with k=0. I need to double-check the reason with Christian Gasser directly, but I surmise that the authors might have used different material properties for that case than those reported in the paper.

In any case, I hope that this response addresses your concerns. For values of k≠0, the agreement between ABAQUS and FEBioStudio is very good. As the authors of the GOH paper indicate, the case k=0 was used for testing purposes but is not a realistic description of the response of the arterial wall to loading.

Best,

Gerard

[athpa]

The agreement between FEBio and ABAQUS for the case k=0.226 is excellent. The relatively small discrepancy in the response for k=0 is, in my opinion, a problem with ABAQUS. The reason I say this is because the case k=0 is equivalent to having two single fiber families. In FEBio, this can be reproduced using an uncoupled solid mixture with a Mooney-Rivlin ground matrix and two fiber-exp-pow-uncoupled fibers ( [ATTACH]n22729[/ATTACH] [ATTACH]n22728[/ATTACH] ). The results for the mixture analysis are identical to the results for the HGO model with k=0, so that gives me confidence that our implementation of HGO is correct.

To enforce strict incompressibility, I used the three-field method in the attached FEBioStudio files. So, if you plot the relative volume J as a function of deformation, you will find that J remains equal to 1 (within at least three decimal places) in the FEBio analyses. I don't know if ABAQUS enforces incompressibility that well, so it wouldn't hurt to check that on your end if you can.

As a final note, the ABAQUS solution does not agree with the corresponding figure published in the Gasser, Ogden, Holzapfel (2006) paper any better than the agreement between ABAQUS and FEBio shown above. Here is a plot of the ABAQUS results superposed over the GOH paper:
ABAQUS-GOH-paper.png

The discrepancy is greatest for the axial speciment with k=0. I need to double-check the reason with Christian Gasser directly, but I surmise that the authors might have used different material properties for that case than those reported in the paper.

In any case, I hope that this response addresses your concerns. For values of k≠0, the agreement between ABAQUS and FEBioStudio is very good. As the authors of the GOH paper indicate, the case k=0 was used for testing purposes but is not a realistic description of the response of the arterial wall to loading.

Best,

Gerard

Hi Gerard,

Now everything is clear. I appreciate your reply.

Regarding the incompressibility and the three-field method , I checked my .feb files and it is missing and I should add it. However, I couldn't find the three-field-solid in the manual, could you share a link so I will have a look?

In Abaqus, to impose strict icompressibility we use the mixed hybrid elements (-H).



Thank you for your time.


Best regards,
Athanasios

[ateshian]

Hi Athanasios,

You can find it described in the User's Manual here.

Best,

Gerard

[athpa]

Thank you Gerard.

Best,

Athanasios

[ateshian]

I have communicated with Prof. Christian Gasser and he suggested that I should double the matrix stiffness parameter c used in FEBio to better reproduce the results of ABAQUS and his paper. When I did that (doubled the value of c used in the files that I posted above), the agreement between FEBio and ABAQUS became nearly perfect for all cases, as shown below:
FEBio-ABAQUS-HGO-model_cx2.png
Best,
Gerard

[athpa]

I have communicated with Prof. Christian Gasser and he suggested that I should double the matrix stiffness parameter c used in FEBio to better reproduce the results of ABAQUS and his paper. When I did that (doubled the value of c used in the files that I posted above), the agreement between FEBio and ABAQUS became nearly perfect for all cases, as shown below:
FEBio-ABAQUS-HGO-model_cx2.png
Best,
Gerard

Hello Gerard,

Amazing ! Thank you for update.

I feel relieved cause we used this publication as a validation tool for our project here in the lab.

I was wondering if FEBio is suitable to predict double anisotropy. What I mean is that we have fiber reinforced matrix(protein) composites that behave anisotropically. The tricky part so far is that the matrix is also anisotropic in both (radial and circumferential) directions.

Best regards,
athpa

[ateshian]

Hi Athpa,

Yes, you can choose to make the ground matrix (e.g.) transversely isotropic, using one of the various trans iso materials (all of which are uncoupled) or one of the orthotropic materials (some uncoupled, some unconstrained -- i.e., compressible), reduced to the special case of transverse isotropy.

You can do this by including the Holzapfel model in an elastic mixture, setting its ground matrix modulus to zero (c=0) and adding one of the above trans iso or ortho materials to the mixture.

Best,

Gerard