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

Collapse
X
 
  • Time
  • Show
Clear All
new posts
  • zigadonik
    Junior Member
    • Nov 2021
    • 16

    #1

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

    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
    Attached Files
  • ateshian
    Developer
    • Dec 2007
    • 2034

    #2
    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

    Comment

    • zigadonik
      Junior Member
      • Nov 2021
      • 16

      #3
      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

      Comment

      • ateshian
        Developer
        • Dec 2007
        • 2034

        #4
        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

        Comment

        • athpa
          Member
          • Aug 2022
          • 36

          #5
          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
          Attached Files

          Comment

          • ateshian
            Developer
            • Dec 2007
            • 2034

            #6
            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

            Comment

            • athpa
              Member
              • Aug 2022
              • 36

              #7
              Originally posted by ateshian View Post
              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 ​
              Attached Files

              Comment

              • ateshian
                Developer
                • Dec 2007
                • 2034

                #8
                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.

                Comment

                • ateshian
                  Developer
                  • Dec 2007
                  • 2034

                  #9
                  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
                  ​​​

                  Comment

                  • athpa
                    Member
                    • Aug 2022
                    • 36

                    #10
                    Originally posted by ateshian View Post
                    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

                    Comment

                    • ateshian
                      Developer
                      • Dec 2007
                      • 2034

                      #11
                      Hi Athanasios,

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

                      Best,

                      Gerard

                      Comment

                      • athpa
                        Member
                        • Aug 2022
                        • 36

                        #12
                        Thank you Gerard.

                        Best,

                        Athanasios

                        Comment

                        • ateshian
                          Developer
                          • Dec 2007
                          • 2034

                          #13
                          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

                          Comment

                          • athpa
                            Member
                            • Aug 2022
                            • 36

                            #14
                            Originally posted by ateshian View Post
                            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

                            Comment

                            • ateshian
                              Developer
                              • Dec 2007
                              • 2034

                              #15
                              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

                              Comment

                              Working...
                              X