×
INTELLIGENT WORK FORUMS
FOR ENGINEERING PROFESSIONALS

Log In

Come Join Us!

Are you an
Engineering professional?
Join Eng-Tips Forums!
  • Talk With Other Members
  • Be Notified Of Responses
    To Your Posts
  • Keyword Search
  • One-Click Access To Your
    Favorite Forums
  • Automated Signatures
    On Your Posts
  • Best Of All, It's Free!
  • Students Click Here

*Eng-Tips's functionality depends on members receiving e-mail. By joining you are opting in to receive e-mail.

Posting Guidelines

Promoting, selling, recruiting, coursework and thesis posting is forbidden.

Students Click Here

Jobs

2D Orthotropic Hyperelastic Fung Model

2D Orthotropic Hyperelastic Fung Model

2D Orthotropic Hyperelastic Fung Model

(OP)
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!

RE: 2D Orthotropic Hyperelastic Fung Model

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.

RE: 2D Orthotropic Hyperelastic Fung Model

(OP)
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.

RE: 2D Orthotropic Hyperelastic Fung Model

(OP)
@ 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!

RE: 2D Orthotropic Hyperelastic Fung Model

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

RE: 2D Orthotropic Hyperelastic Fung Model

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?

RE: 2D Orthotropic Hyperelastic Fung Model

(OP)
@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:



RE: 2D Orthotropic Hyperelastic Fung Model

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

RE: 2D Orthotropic Hyperelastic Fung Model

(OP)
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:

RE: 2D Orthotropic Hyperelastic Fung Model

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.

RE: 2D Orthotropic Hyperelastic Fung Model

Quote (Nucleophobe)

"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.

Quote (Nucleophobe)

"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

RE: 2D Orthotropic Hyperelastic Fung Model

(OP)
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.

RE: 2D Orthotropic Hyperelastic Fung Model

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.

RE: 2D Orthotropic Hyperelastic Fung Model

(OP)

Quote (sdebock)

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/jAfnfj1...

RE: 2D Orthotropic Hyperelastic Fung Model

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.

RE: 2D Orthotropic Hyperelastic Fung Model

(OP)
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:

Quote (Abaqus)


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

RE: 2D Orthotropic Hyperelastic Fung Model

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

Quote (Abaqus)

***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?

RE: 2D Orthotropic Hyperelastic Fung Model

S is not Cauchy
In abaqus it is either:
2nd Piola Kirchoff stress ( = J*F^(-1)*sigma*F^(-1) )
or:
The deviatoric stress (S = sigma + p * I)

In this case, I'm pretty sure it's the deviatoric stress.
So to get Cauchy you do
sigma_ii = -p (=1/3 sig11+sig22+sig33) + lambda_i* dW/dlambda_i

Meaning (I guess) your analytical solutions were 1/3th to high?

And yeah, Shell elements can have incompressiblity (i never checked before). Because they are reduced integration, normally you shouldn't get volumetric locking either.

RE: 2D Orthotropic Hyperelastic Fung Model

(OP)
Not yet. I'm reading Holzapfel's book over Christmas break so I understand the fundamentals better.

I found other sources saying that S is Cauchy, but I need to check the documentation.

RE: 2D Orthotropic Hyperelastic Fung Model

Correct, ABAQUS uses Cauchy stress.

By the way, I finally got the chance to look at your INPs and there is one major concern with the boundary conditions: You are not constraining the rotational degrees of freedom at all nodes. As you know, nodes on shell elements (well, structural elements in general) have both translational and rotational degrees of freedom. Without constraining those DOFs, you are not enforcing a uniaxial stress state, which makes it impossible to do a reasonable comparison with theoretical predictions.

Another concern: kappa (the fiber dispersion parameter) seems to be quite high. Is this value of kappa relevant to the histological/fibrous nature of the tissue? Keep in mind that you are already assigning two directions of anisotropy (fiber directions; you are not explicitly assigning discrete fibers, just directions and degrees of anisotropy) which are at some angle with each other. Make sure you have your angles defined correctly.

Note that ABAQUS will treat "Axis 1" as the primary/longitudinal/along-fiber loading direction whereas other directions will be orthogonal to it (unless you provide some additional rotation or a spatially varying orientation using discrete field).

Practical notes about k1 and k2 (check this in your MATHEMATICA notebook):

k1 - (stress-like parameter)

- smaller k1 will create a longer toe region (imagine collagen fibrils being too wrinkled to resist any load)
- higher k1 shortens the toe region

k2 -

- higher k2 creates a sharper increase in stress (i.e., higher number of collagen fibrils get recruited to resist applied load)
- smaller k2 creates a smoother increase in stiffness (so, a very small k2 will increase the toe region).

Hope this brings some joy winky smile

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

RE: 2D Orthotropic Hyperelastic Fung Model

Maybe we got mixed up. I thought you were asking about the S used in documentation.
Abaqus indeed uses Cauchy, but in documentation manual, sigma symbol is used for Cauchy, and S for 2PK or dev stress.
In CAE, S is Cauchy.

RE: 2D Orthotropic Hyperelastic Fung Model

(OP)
Thank you both! Yes, I meant to ask what 'S' in Abaqus CAE represents. Thank you for confirming sdebock.

I'll try constraining the rotational degrees of freedom as suggested by IceBreakerSours. I didn't think that a uniaxial displacement would have any influence on the node rotations.

RE: 2D Orthotropic Hyperelastic Fung Model

(OP)
Hey again,

I added the constraints for the rotational degrees of freedom where they were missing as IceBreakerSours suggested. The 'kappa' I am using is consistent with the example I found in the documentation. The mathematical model also seems to work okay with this value.

I'm still getting the same results. Do you know if I am defining my invariants 'I1' and 'I4' correctly? I am using equations that I found in the Holzapfel paper Hyperelastic modelling of arterial layers with distributed collagen fibre orientations (2006). This is very close to the form cited in the Abaqus documentation.

RE: 2D Orthotropic Hyperelastic Fung Model

I ran your input file.

your lambda_t goes till 1.3.
from output, we see that lambda_x is 0.654688, which is not equal to the sqrt(1/lambda_t) it should be from your formula.
What's your thought on this?

RE: 2D Orthotropic Hyperelastic Fung Model

(OP)
@sdebock,

Thank you for trying the input file! By lambda_t, do you mean the theta direction? I stretch the 10x10mm element 1.5mm in each direction (one per file), which should yield a lambda_t of 1.15. Can you clarify why you get a lambda_t of 1.3?

How did you get lambda_x from the output? Do you just query the width of the element in the x direction in the deformed state and calculate? I would like to check this again. Thank you for the pointer.

Does this suggest that Abaqus is not treating the material as fully incompressible?

By the way, I updated the input files to simplify the problem. Before I was both compressing and stretching the element. Now I only do the compression step. The sets are also named.

RE: 2D Orthotropic Hyperelastic Fung Model

(OP)
UPDATE:

Yes, I believe this is the problem!

I made the assumption that in uniaxial tension in the lambda_1 direction, lambda_2 = lambda_3. However, since the material is anisotropic, my intuition suggests this to be incorrect. The constitutive equation for incompressibilitiy should still hold however:
labmda_1*lambda_2*lambda_3 == 1
OR
labmda_t*lambda_r*lambda_z == 1

The problem is, how do I relate the lambda's? In order to use the starting SEDF, I need to have:
1) For uniaxial tension in the theta direction: lambda_z as a function of lambda_t
2) For uniaxial tension in the z direction: lambda_t as a function of lambda_z

Yikes.

RE: 2D Orthotropic Hyperelastic Fung Model

Hooray, I have solved your problem.
Incompressibility still holds.
So, if you hold lambda_Z to be unknown, you have a problem.
Luckily, you have a spare equation hanging around doing nothing, like a chump.
If you pull in theta direction, you solve for sigma_theta.
The trick now is to solve sigma_z = "long equation with only lambda_theta and lambda_z unknown" == 0 !
And use this lambda_Z = x*lambda_theta in sigma_theta = "...."
My math program was not able to solve this equation by pressing solve, but I've confirmed it for some numbers of lambda_theta & abaqus.
But surely with mathematica or whatever you should get an analytical answer & a nice curve.

RE: 2D Orthotropic Hyperelastic Fung Model

(OP)
Excellent! I never thought to use sigma_2 == 0 when solving for sigma_1 as a function of lambda_1! Yes, Mathematica should have no problem with this.

I can't wait to try this. I owe you both, you have been incredibly patient and helpful!

RE: 2D Orthotropic Hyperelastic Fung Model

Also, why are you using the cylindrical orientation for a square uniaxial test strip?
I confirmed it for C3D8H elements, since I have most experience with those.

RE: 2D Orthotropic Hyperelastic Fung Model

(OP)
I'm using cylindrical because the final mesh will be round (see images above). Using cylindrical coordinates, I don't have to define a local material orientation for each element; I simply assign the coordinate system and move forward.

I'm doing it in the validation just to make sure it works as I expect it too. You have a point though; I should probably try it in Cartesian first.

RE: 2D Orthotropic Hyperelastic Fung Model

In this thread I state:

""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 brackets! ""

This is actually not true (my bad), but depends on the value of stretch, gamma & kappa! So you should always check it!

RE: 2D Orthotropic Hyperelastic Fung Model

Nucleophobe:

Can you upload your latest INP? Also, let us know what you have tried recently and what corresponding issues you are facing with the current model. It shouldn't have taken this long to get the model to work.

sdebock:

If you weren't referring to the Mathematica notebook, then while the amount of stretch, fiber orientation, and fiber distribution about a mean orientation all determine what the stress is going to be (besides the strain state), the term with Macaulay bracket is to allow for fibers to not take any compression.

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

RE: 2D Orthotropic Hyperelastic Fung Model

(OP)
Eureka! I just finished implemented the changes to the Mathematica code and checking against Abaqus. Using the 'sigma-transverse == 0' equations, the system is harder to solve, so I have to work out a new method for fitting the parameters. I can deal with that however.

@IceBreakerSours:
I uploaded my current files. Here's the link (same as before):
https://www.dropbox.com/sh/8imp0znuspu2tzo/jAfnfj1VJ3

The PDF is also updated.

@sdebock:
I took your advice and verified E_alpha to be positive for all lambda values of interest.

I hope this helps others too. Let me know if you have questions.

RE: 2D Orthotropic Hyperelastic Fung Model

(OP)
I'm done, thanks again!

Note: my directions are backwards in the analytical model; I'm checking to see why, but not a big problem.

Red Flag This Post

Please let us know here why this post is inappropriate. Reasons such as off-topic, duplicates, flames, illegal, vulgar, or students posting their homework.

Red Flag Submitted

Thank you for helping keep Eng-Tips Forums free from inappropriate posts.
The Eng-Tips staff will check this out and take appropriate action.

Reply To This Thread

Posting in the Eng-Tips forums is a member-only feature.

Click Here to join Eng-Tips and talk with other members!


Resources