convergence and fibre orientation

[vittoria]

the first one!

I have downloaded FEBio and I have started using it.

I think it has amazing capabilities so I have started using it for my PhD project.

I conducted a uniaxial and a biaxial test using a transverse isotropic Mooney-Rivlin material model.

I used my constants and also the ones found on Jeff Weiss’ thesis. What I experienced is that both simulation were not converging even if I was using for small forces (0.01 N or even 0.00001N). They would converge only for stretches lower than 1.1.

FEBio should be a nonlinear finite element software package but with those low stretches and low forces I can not study the non-linear behavior I am interested in.

What can be the problem? Where I am wrong?

I do not think are my constants because also with the ones from the thesis it does not work. Moreover, I think the constraints are ok.

Furthermore, I have a doubt about the fiber local orientation using nodes. If I want them to be oriented on z, what do I have to write: 1, 5? Can I orient the fibers using all the nodes of a hexahedral element or there is some restriction?

Thanks for all the support you will give me

v

[maas]

RE: convergence and fibre orientation

Hi Vittoria,

After looking at your files I believe that the problem does lie with your material properties. For example, in the biaxial model you report a value for C4 of 7.69e+007. However, as far as I can tell from the literature this value is usually less than 100. Here are some usefull links with material parameters: Quapp1997 and gardiner2003. For example, when using the first set of material parameters from table 1 in the second article your problem runs without any issues.

Some other remarks about your files. In the uniaxial model your loadcurve extends the total timerange. In the biaxial model the y-force and z-force are of by a factor of 1e6. I have attached my modified files. These files should run without problems.

As far as your other fiber orientation question is concerned, you can use all eight nodes to specify the orientation of the fiber. However, if you want to orient all vectors along the z-axis it might be easier to use the "vector" option. For instance, the following line will orient all fibers along the z-direction:

<fiber type="vector">0,0,1</fiber>

This is easier since it does not require explicit knowledge of the local element node numbering. I hope this helps.

Steve.

[vittoria]

Hi Steve,

thanks for the exhaustive and quick answer.
Now, following your suggestions, I am trying to run my uniaxial test with the data from the first paper (table 4, page 17, first row).
What happens is that it converges, but only for small loads. The highest stretch that I can get is 1.003, that is low.
This is low also when compared with the uncrimped stretch, that is of 1.035for this dataset.
Can I ask you if there is any reason for this? Is there any parameter in the solution solver that I can change?

Thanks again

Vittoria

[maas]

Hi Vittoria,

Can you please attach the file you are having trouble with.

Thanks,

Steve.

[vittoria]

Sorry for the delay, but I was abroad and had not the files with me.

Anyway here are the files. One has a coarse mesh and the second a refined one.
They both do not converge, and I think it is strange as the maximum stretch reached is in both cases around 1.003, while the lambda max in the material model is 1.035.
I think I should be able to reach the lambda max of the material model.

I run the files increasing the load. I am also attaching an image of the uniaxial test for the highest converging load, where you can see the level of strain reached.

Thanks

Vittoria

[maas]

Hi Vittoria,

I believe the problem is that your bulk modulus (k) is too small. If I set it in both cases to 100MPa I can apply much larger loads. (~ x1000). Also, keep in mind that using negative values for C2 may in some cases lead to instabilities. It is for this reason that you'll notice that in the second paper I reference above they've used C2=0 in all their analyses.

Another thing I want to point out is that you apply the same force on all top nodes. It may sound strange, but this will not generate a uniform force at the top as you might expect. The reason for this odd behavior lies in the details of the FE method. One way of understanding it is that nodes that are attached to more elements will appear stiffer and move less. You'll see what I mean once you start applying larger forces. In the next release of FEBio we'll have a way around this: You can attach the nodes to a rigid body and apply the force on the rigid body. Let me know if you have any other questions.

Cheers,

Steve.

[vittoria]

Thanks for all, now it is working.

I noticed the strange behaviour applying the same load to all the nodes...but how can I know how much the nodal load has to differ from node to node? Is there a relationship based on the number of elements the node is attached to? (say x4 if the elements are for or something like that)

For now I am using a pressure, but it converges only for small strains. I am attaching the file.

Thanks again

Vittoria

[maas]

Hi Vittoria,

After running your problem on a different platform (using a different solver) and comparing the results with a different program I am begginning to suspect a bug in one of the stiffness routines. The good news is, is if it converges the results are identical to this other program. The bad news is that it should run for a lot bigger pressures. I'll look into it more and let you know what I find.

Thanks,

Steve.

[maas]

Hi Vittoria,

I believe I have found the source of the problem. It appears that the pressure stiffness term is creating instabilities in your analysis. The pressure stiffness term is added to the global stiffness matrix for problems using prescribed normal pressures (like in your case) and is supposed to help the problem to converge quicker. However, quit the opposite is going on in your case (which might have something to do with the fact that the pressure boundary is free, although I'm not sure). But here's what I'll do. I'll add an option to FEBio that allows you to turn this pressure stiffness on or off. I'll default it to off since it appears most problems are running better without this pressure stiffness. Also, that way you won't have to make any changes in your input files.

Cheers,

Steve.

[Adam Baker]

possibly related

I don't know whether this is related to the above problem or not. Previously, when I ran the attached code, I got a negative Jacobian error (and WinXp, Linux32, and Altix), but with 1.1 it completes successfully. Nothing required since it works now, but I thought I would post the code in case it was of any assistance.

Adam

[vittoria]

Hi,

Thanks again for all the support you previously gave me.

Now I overcame all the past problem and I started using the transversely isotropic hyperelastic constitutive model for different fiber angles (using forces and not pressure!).

I have been changing the fiber angle using the command you suggested me (<fiber type="vector">0,0,1</fiber>) and changing the values according to the sin and cos of the fiber angle.

What I obtain is that for fiber angles bigger than 30 degrees to the pulling direction all the strain-stress curves flattens.
To make it clearer there is no difference between a curve obtained for fibers perpendicular to the pulling direction and the one obtained for fibers at 40 degrees to the pulling direction.

This seems strange to me as I thought that the fibers at intermediate angles were still giving a contribution to the stress strain curve.

I am attaching the code written for a uniaxial tensile test with fibers oriented at 30 degrees to the pulling direction.

I hope I have exposed clearly what my problem is.

best regards

Vittoria

[maas]

Hi Vittoria,

So, I've ran your model for the case 0 and 90 degrees and compared them with the analytical resulsts. I have attached my findings and as you can see they agree quite well. The 30 degree case did not correspond to my analytical result but after looking at your result I realized that this is probably because your mesh shears for this particular case and has a non-uniform deformation. My analytical result on the other hand assumes perfect uni-axial tension. To compare to the analytical result we'd somehow have to setup boundary conditions to enforce a uniform deformation in the mesh, which now that I think about it, is hard to do. I'll give it some more thought and get back to you once I've come up with a solution. Hopefully, the good correspondence between the analytical results for the 0 and 90 degree cases will give you some better confidence in the FEBio results.

Cheers,

Steve.

[vittoria]

thanks a lot!

anyway I was confident with FEBio results!

[vittoria]

Hi Steve,

I have been trying to change the constraints, but there are still no improvements.
I was wondering if you could send me also a plot of the curve obtained for 30 degrees using your analytical results, so that I can compare it with my results.

Thanks again

Vittoria

[shore]

Steve,

Is there any way to use PreView to determine the local element node numbering and/or what is it by default? I built my geometry in PreView (a bunch of cylinders built from the bottom up) and I'm trying to apply the transversely isotropic Mooney-Rivlin model with fibers running around the circumference (its an intervertertebral disc annulus). I assume I can use local n1 = 1 and n2 = 2 for the fiber orientation, I just want to verify that my nodes are numbered right side up. Thanks.

Spencer Shore