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!

 

[burningDiesel]

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.

 

[IceBreakerSours]

Did you manage to solve your problem(s)?

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

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

 

[Nucleophobe]

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.

 

[IceBreakerSours]

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 ;-)

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

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

 

[burningDiesel]

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.

 

[Nucleophobe]

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.

 

[Nucleophobe]

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.

 

[IceBreakerSours]

Some general comments:

a) ABAQUS/CAE provides the option of defining the constants for the HGO model directly, there's no need to worry about invariants.
b) Make sure the orientation of 1-axis is what it is supposed to be.
c) Make sure the angle you define is correct.
d) Check your units.

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

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

 

[burningDiesel]

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?

 

[Nucleophobe]

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

 

[Nucleophobe]

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.

 

[burningDiesel]

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.

 

[Nucleophobe]

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!

 

[burningDiesel]

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.

 

[Nucleophobe]

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.

 

[burningDiesel]

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!

 

[IceBreakerSours]

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

 

[Nucleophobe]

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):
[URL unfurl="true"]https://www.dropbox.com/sh/8imp0znuspu2tzo/jAfnfj1VJ3[/url]

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.

 

[IceBreakerSours]

Do you still need help with implementing the model or are you done?

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

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

 

[Nucleophobe]

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.