2D Orthotropic Hyperelastic Fung Model

[Nucleophobe]

I am working on modelling the response of an orthotropic hyperelastic material (vein) to a given load. Because the local thickness of the material is unknown, shell elements with a constant thickness are used. I am able to get the simulation to run. However, I had some questions about the material properties.

The orthotropic material in Abaqus requires 11 material constants: b1111, b2222, b3333, b1122, b1133, b2233, b1212, b1313, b2323, c, and D. The parameter 'c' defines the linear response of the material, while 'D' determines the response to a thermal load and may be set to 0.

Since this is a 2D shell model, I thought all terms involving the '3' direction would drop out. Abaqus requires that b3333, b1133, b2233, b1313, and b2323 are non-zero however. I have been setting them to arbitrary values and observing the results. The solutions seem to be quite similar with widely varying properties in the '3' direction; however, the solutions do not converge when these constants differ from my b1111 and b2222 constants by an order of magnitude or more.

Needless to say I'm a little lost concerning the theory behind this material model. With a 2D orthotropic material, shouldn't I only need the modulus in the '1', '2', and '12' directions? That would leave b1111, b2222, b1122, and c.

Any comments / recommended references would be greatly appreciated. Thank you!

 

[MechIrl]

I'm not so sure about the Abaqus implementation of shells. However, I think the 'out-of-plane' properties are used by Abaqus to determine the shell's 'transverse shear stiffness'. Also 'D' (as it is in the Abaqus manual) is not the response to thermal load in the Fung model. It is a measure of the bulk modulus of the material and is only zero if the material is incompressible.

 

[Nucleophobe]

Thanks Mechlrl, that helps. I saw that 'D' was temperature dependent in the documentation and didn't interpret that correctly. Incompressibility is often assumed for veins since they are mostly water; now I can justify using 0 correctly.

The root of the issue is that biaxial test data for this particular vein (and most veins in general) is lacking. I am trying demonstrate how the material could be modelled by fitting two uniaxial stress-strain curves to Fung-orthotropic material constants. I guess I will need to substitute *reasonable* parameters for the out-of-plane behavior in order to define the 'transverse shear stiffness'. I'm having trouble finding any literature where a Fung-orthotropic model is used actually (most use Holzapfel, but the Abaqus implementation of Holzapfel yields a hyperelastic response only in one direction and I was hoping to avoid issues with fiber orientation).

On a related note, do you know how to represent the material constants b_ijkl in a stiffness matrix? I sometimes get the error:
"The stiffness matrix defined for the fung-orthotropic hyperelastic material is not positive definite"
but I'm not sure how this is computed.

 

[IceBreakerSours]

I am not sure why you are not keen on using HGO but if you allow two directions of anisotropy and have them at appropriate orientations, then HGO material model should give you a reasonable fit against the experimental data.

http://www.eng-tips.com/faqs.cfm?fid=376

http://www.eng-tips.com/faqs.cfm?fid=1083

 

[Nucleophobe]

@ IceBreakerSours:

The reason I moved away from the HGO model is that I had trouble getting two material orientations for the fibers using CAE.
When you say "allow two directions of anisotropy", do you mean specifying "Number of local directions" to be two? Would I then specify two sets of material coordinates, one for each of the fiber directions?

The model has a couple bifurcations, so specifying material coordinates is a little tricky. I have been doing so using a cylindrical datum coordinate system with the origin roughly in the center of the vein and the R-theta plane cutting perpendicular to the vein. I then assign an orientation using this coordinate system, which automatically assigns a '1' and '2' direction to the elements in a selected region as desired.

Can I set the fiber directions to be based on these coordinates, but rotated +-30degrees based on the unit normal? That would be great!

 

[IceBreakerSours]

You can not assign directions of anisotropy (i.e., fiber orientations) in CAE; it must be done in the input file directly. Check the input files in the Examples, Verification and other manuals for details. In particular, check the parameteric input file(s) that assign an angle between two fiber families.

What do you mean by "bifurcations"? Yes, you can assign two directions of anisotropy (there are no "discrete fibers", just orientations of anisotropy) in a straightforward manner. Check the Analysis Manual along with the associated input files.

http://www.eng-tips.com/faqs.cfm?fid=376

http://www.eng-tips.com/faqs.cfm?fid=1083

 

[IceBreakerSours]

Edit: It is much more convenient to assign directions of anisotropy in the input file directly than CAE.

http://www.eng-tips.com/faqs.cfm?fid=376

http://www.eng-tips.com/faqs.cfm?fid=1083

 

[MechIrl]

I'm not very familiar with the Fung model so can't really help much regarding the stiffness matrix. If you can get your hands on Abaqus 6.12 the HGO model is available in CAE, maybe it will make life easier since the 'out of plane' behavior will be automatically given by the isotropic hyperelastic part of the model.

Dealing with bifurcations with the HGO model is tricky but like IceBreaker says there are a couple of methods available to assign directions when you use the INP fie. This includes something along the lines of what you were looking for in terms of allowing additional rotations about a given co-ordinate system. The OFFSET-TO-NODES option might also be useful as (in 3-D elements at least) it allows the element's co-ordinate system to be defined in terms of its node's positions.

I also imagine you may face similar difficulties with the Fung model in terms of assigning directions since the material is othrotropic?

 

[Nucleophobe]

@IceBreakerSours
By "bifurcations", I mean that the vein branches. I should actually call them anastomoses I guess (joining of two smaller veins to form a larger vein). I've been defining a local cylindrical coordinate system for each branch up to now, which has worked out okay. I'll see how it goes trying this approach from the input file.

@Mechlrl
I was using the Fung model because I have uniaxial test data in the circumferential and longitudinal directions. It was fairly easy to fit parameters to the '1' (longitudinal) and '2' (circumferential) directions. I then assigned material orientation in CAE using a cylindrical coordinate system as explained earlier.

It sounds like doing this will be much more difficult with the HGO model. Are there any examples of the HGO model with two directions of anisotropic hyperelasticity? Whenever I use it, I get a nonlinear stress-strain curve in the '2' direction but not in the '1'. Does this mean I need to define two 'local directions', which would represent the collagen fibers wrapping left-hand and right-hand around the vein (like left-hand and right-hand screws)? I will then have to mess with the angle until I get the right stress-strain response in the two directions... I wish my composites book hadn't gotten lost in the mail when I moved!

Here are some screenshots. You can see that each element has a local coordinate system as desired, and it only takes a minute or so to setup:

 

[Nucleophobe]

Edit: I made this simplified model to explain. The actual model includes multiple branches and is not as uniform.

 

[Nucleophobe]

Edit 2: I'm looking in the 6.12 documentation for HGO now to see how the fiber orientations are specified.

Here is the desired behavior of the material model. You can see that both curves are highly non-linear, with the longitudinal direction about twice as distensible as the circumferential direction:

 

[MechIrl]

Yes you will need two fibre directions for the HGO model for your curves. The paper 'Determination of layer-specific mechanical properties of human coronary arteries with nonatherosclerotic intimal thickening and related constitutive modeling' describes how to fit a version of the HGO model (a little different to the one in Abaqus) to very similar data to what you have.

Really its not such a big deal to assign fiber directions once the local coordinate system is known. If you already have this system for you elements in the bifurcation you've already done the hard part.

 

[IceBreakerSours]

"Are there any examples of the HGO model with two directions of anisotropic hyperelasticity?"

There are many such examples in the manual. If I remember correctly, there are a few examples that verify results from one of Gasser's papers on arteries. In those examples, ABAQUS employs (as in the paper) two fiber families at an angle with each other.

"Whenever I use it, I get a nonlinear stress-strain curve in the '2' direction but not in the '1'."

The angle between fiber families must not be set correctly. If you read the Analysis User's Manual carefully, you'll find that setting the angle is a piece of cake. If you have the code/script for curve fitting ready, as MechIrl noted, the hard part is done.

http://www.eng-tips.com/faqs.cfm?fid=376

http://www.eng-tips.com/faqs.cfm?fid=1083

 

[Nucleophobe]

Thank you all again, I have made progress using the HGO model. I'm still having trouble fitting parameters to my data. The form of the strain-energy density function seems a little weird to me.

If I understand correctly, the Cauchy stress 'sigma' in the 11 direction would be:
sigma_11 = lambda_1 * partial(U)/partial(lambda_1)

where 'U' is the strain energy density function as defined in the Abaqus documentation (neglecting the compressibility term):
U = C_10(I_1-3) + k_1 / (2 k_2) * Sum[ exp [ k_2 * <E_alpha>^2 ] - 1 ]

However, when taking this derivative, the <E_alpha> term must be dealt with. The documentation mentions that <-> stands for the Macauley bracket, which performs the operation:
<x> = 1/2(Abs[x]+x)

I'm not sure how to differentiate this appropriately.

Do you have any tips regarding how to get the uni-axial stresses sigma_1 and sigma_2 from the strain-energy density function 'U' used by Abaqus? The Yahoo! group has a discussion on this, but it sort of dead ends.

 

[burningDiesel]

The <x> term just means that the fibers don't act under compression. Since u are fitting uniaxial tensile data, you can just ignore the bracket!
If you want to continue, the derivative of the ramp function is the 'Heaviside step function' I guess.

 

[Nucleophobe]

The <x> term just means that the fibers don't act under compression. Since u are fitting uniaxial tensile data, you can just ignore the bracket!
If you want to continue, the derivative of the ramp function is the 'Heaviside step function' I guess.

Ah yes, E_alpha is positive for all lambda except for extreme cases of gamma, so that should be okay.

In that case I must be doing something else wrong; I'm verifying against Abaqus using a single 10x10mm shell element in tension. I displace one side of the element by 1.5mm, yielding a stretch of 1.15. I do the same thing again in the opposite direction and plot them to compare with my analytical model. The stress-strain curves resulting from Abaqus are more shallow than those predicted by my code.

Here are the files if you're willing to look at them. Holz-D1.inp is the circumferential direction, and Holz-D2.inp the longitudinal. The code is in Mathematica, but there's a PDF as well. If nothing else, might help someone else:
DropBox

https://www.dropbox.com/sh/8imp0znuspu2tzo/jAfnfj1VJ3

 

[burningDiesel]

disclaimer: I didn't check your maths, only your input files.
Assuming your analytics are correct:

Why do you check it on a shell element?
First do a uniaxial tensile test on a single C3D8H element, the H stands for incompressibility. Right now, you are using D=0, but Abaqus can not handle a poissons = 0.5 unless you use an incompressible element. Probably in your .dat file, you are getting the warning: "Abaqus set the value of D to xxx.xxx", which corresponds with a poissons of 0.475 (I think).

So either put the D that you get in that warning in your analytical model, or use the incompressible element.

If there is only a small difference between your abaqus & analytical curve, this might be it.

 

[Nucleophobe]

I checked both the .msg and .dat files and couldn't find the suggested warning. Also, in the Abaqus documentation example, they suggest that D=0 is acceptable:

Holzapfel-Gasser-Ogden energy function coefficients:
C10 = 3.82 kPa
k1 = 996.6 kPa
k2 = 524.6
kappa = 0.226
D = 0 ( = 1× 10–6 for compressible case)

Anyways, small changes to the compressibility do yield large changes in material response. Have you seen this paper? Its very recent, and claims that there are serious problems with HGO:
"Deficiencies in numerical models of anisotropic nonlinearly elastic materials" (Annaidh, Sept. 2012)
Link

 

[IceBreakerSours]

It may be that incompressibility constraint is preventing the nodes to move as freely as they would otherwise; the problem, in its current form, may be "too stiff". Instead of using a single shell element, try using multiple elements.

http://www.eng-tips.com/faqs.cfm?fid=376

http://www.eng-tips.com/faqs.cfm?fid=1083

 

[Nucleophobe]

I figured out which elements may be used with incompressibility by trying a C3D8 element and getting this error:

***ERROR: INCOMPRESSIBLE ANISOTROPIC HYPERELASTIC MATERIALS CAN ONLY BE USED WITH HYBRID, PLANE STRESS, OR SHELL ELEMENTS

Looks like shell elements are okay. The final geometry we will be using is a 3D surface with uniform thickness, so shell elements are the most convenient. That's why I've been trying to use shells to validate the analytical model.

@IceBreakerSours:
I'll try what you suggested, though the Abaqus results are actually "too compliant". More compliant stress-strain results would make sense if Abaqus were assuming some incompressibility, but apparently it isn't.

I suspect something is wrong with the math... Or my interpretation of nomenclature. 'S' in Abaqus is Cauchy stress, right?