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!
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
RE: 2D Orthotropic Hyperelastic Fung Model
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
http://www.eng-tips.com/faqs.cfm?fid=376
http://www.eng-tips.com/faqs.cfm?fid=1083
RE: 2D Orthotropic Hyperelastic Fung Model
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
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
http://www.eng-tips.com/faqs.cfm?fid=376
http://www.eng-tips.com/faqs.cfm?fid=1083
RE: 2D Orthotropic Hyperelastic Fung 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
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
RE: 2D Orthotropic Hyperelastic Fung Model
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
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
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.
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
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
If you want to continue, the derivative of the ramp function is the 'Heaviside step function' I guess.
RE: 2D Orthotropic Hyperelastic Fung Model
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
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
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
http://www.eng-tips.com/faqs.cfm?fid=376
http://www.eng-tips.com/faqs.cfm?fid=1083
RE: 2D Orthotropic Hyperelastic Fung Model
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
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
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 found other sources saying that S is Cauchy, but I need to check the documentation.
RE: 2D Orthotropic Hyperelastic Fung Model
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
RE: 2D Orthotropic Hyperelastic Fung Model
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
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
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
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
RE: 2D Orthotropic Hyperelastic Fung Model
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
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
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
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
I can't wait to try this. I owe you both, you have been incredibly patient and helpful!
RE: 2D Orthotropic Hyperelastic Fung Model
I confirmed it for C3D8H elements, since I have most experience with those.
RE: 2D Orthotropic Hyperelastic Fung Model
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
""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
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
@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
http://www.eng-tips.com/faqs.cfm?fid=376
http://www.eng-tips.com/faqs.cfm?fid=1083
RE: 2D Orthotropic Hyperelastic Fung Model
Note: my directions are backwards in the analytical model; I'm checking to see why, but not a big problem.
RE: 2D Orthotropic Hyperelastic Fung Model
By the way, once you resolve the remaining issues, please upload all your files. They may turn out to be beneficial to others who visit this thread.
http://www.eng-tips.com/faqs.cfm?fid=376
http://www.eng-tips.com/faqs.cfm?fid=1083