Axisymmetric Analyses in FEBio

[ateshian]

Note: The contact interface "sliding-tension-compression" becomes available starting with FEBio 1.5.2.

FEBio is a 3D finite element modeler which, as of the writing of this tutorial, does not have 2D axisymmetric elements. This tutorial shows how to simulate an axisymmetric analysis using only one sliver of elements along the circumferential direction.

For illustration, consider a cylindrical elastic solid subjected to a uniform compressive pressure on the top face and a fixed bottom face. This analysis may be performed on the full geometry:
DiskFull.jpg

More commonly, using the x-z and y-z coordinate planes as symmetry planes, the analysis may be performed on a quarter geometry (DiskUniformPressureQuarter.feb):
DiskQuarter.jpg

However, significant computational efficiency would be gained by performing the analysis on a single sliver or wedge portion of the axisymmetric problem (DiskUniformPressureSliver.feb):
DiskSliver.jpg

To perform the latter analysis, the cross-sections on either side of the wedge need to be defined as planes of symmetry.

• The face lying in the x-z plane may be constrained as usual (fixed along y).
• If the geometry includes an edge on the z-axis (as in this example), fix the x and y displacements for that edge.
• To constrain the face lying at an angle to the x-z plane, create a planar surface in the x-z plane and rotate it about the z axis by an angle equal to the wedge angle. Define a rigid body material and attach it to this symmetry plane, then fix all six of its degrees of freedom.
• Define a contact interface of type "sliding-tension-compression" between the symmetry plane and the mating wedge cross-section (single pass, slave surface on the wedge, master surface on the rigid plane, auto-penalty on, penalty=1e4) and make sure to turn on the tension flag (to allow this sliding interface to sustain tension).
• Apply all other relevant boundary conditions on the wedge and run the analysis.

Check the attached file DiskUniformPressureSliver.feb for more details. Even though a contact analysis needs to be performed on the rotated symmetry plane, the analysis runs much faster than the quarter-symmetry model.

Notes:

1. Make sure that the outward normal to the symmetry plane faces opposite to the outward normal of the mating wedge cross-section.
2. In principle, the smaller the wedge angle, the greater the fidelity of the model to true axisymmetry. As the wedge angle is decreased, consider increasing the contact penalty to ensure accurate enforcement of the symmetry requirement.

[Qingen]

Thank you very much, Prof. Ateshian. I believe it is very helpful for researchers needing to perform axisymmetric analyses.

[guerin]

Recommandation about axisymetrical analysis

Hi,
Thank you for this new feature.
I have a question : do you think that modelling the same sliver using two "Rigid wall interfaces" with opposite normal and same slave surface could be a solution? If yes, which one do you recommend as a numerical stability and cost point of view ?
Guerin F.

[ateshian]

Hi Guerin,

A rigid wall interface is like a regular contact interface in the sense that it can only sustain compression, not tension. For axisymmetric problems, it is often the case that there will be tension on the symmetry plane (as in the example of this tutorial). With a rigid wall interface, the deformable material will simply pull away (separate) from the rigid wall in that case. This will violate the symmetry constraint. That's why we needed to introduce a new frictionless contact interface that can sustain tension ("sliding-tension-compression").

Of course, you may be able to define problems (or specific material properties) that produce compression instead of tension on the planes of symmetry. In that case, the rigid wall interface should work. The drawback is that the validity of this approach becomes problem-specific.

Best,

Gerard

[campolee]

Gerard,

We are trying to adapt the sliver technique you describe above (for a solid cylindrical domain) to a biphasic, pressurized hollow spherical domain. We are getting good agreement for the fluid pressure distribution against the analytically solution of Rice and Cleary (1976), but the displacement results are ... funny. Specifically, the y and z displacements are both positive (but unequal in magnitude), but the x displacement is negative, despite the domain residing in the positive octant, as shown in the attached images. Are we specifying something incorrectly in our simulation file (attached)? We have reviewed the boundary conditions we apply against the sample cylindrical feb simulation you provided, but we do not see any obvious error. Any assistance you can provide would be greatly appreciated.

Thanks,
Eamon

[ateshian]

Hi Eamon,

I noticed a couple of things in your file that you should address, though they did not fix the problem you describe:

1) The patch elements you used for the symmetry plane have their surface normals pointing in the wrong direction (away from the spherical wedge). You need to regenerate the patch with the correct orientation of the normals.

2) The penalty factor for the sliding-tension-compression contact interface should be significantly greater than 1, I suggest 1e4.

3) I assume that your intention is to prescribe just a fluid pressure p0 on the inner surface of the hollow sphere. In this case, keeping in mind that the total traction t is equal to t=-p+te, where te is the effective traction, your boundary conditions should be either (p=p0, t=-p0) or (p=p0, te=0). It appears that you applied (p=p0,te=-p0).

Unfortunately, fixing these issues will not overcome the problems you are encountering. It is not just the displacement which seems incorrect in the output, but also the principal stresses don't seem to satisfy spherical symmetry. So I am concerned that there may be something wrong with the mesh of the hollow sphere, though I can't say what exactly. If I were performing this analysis myself, I would generate a mesh with only one element along the azimuth and elevation angles, and multiple elements along the radial direction. The mesh along the radial direction should be biased to produce much finer elements near the inner spherical surface, since the fluid pressure will decay rapidly from there toward the outward spherical surface. For this type of mesh you would need two rigid planes of symmetry, one rotated along the azimuth and the other along the elevation. I would suggest you try that out.

Best,

Gerard

[campolee]

Gerard,

Thank you for your response. We noticed in your original posting the comment about the outward normals, but we do not know how to visualize those in PreView. We tried activating the normals via the checkbox through the display tab of the options menu, but that did not display anything (we are using PreView 1.12.2 for Windows 7).

Otherwise, we will adjust the penalty factor to your recommended value of 1e4 (we did vary that and did not notice any effect) and correct the traction boundary condition. We will also try your suggestion to make a spherical wedge with two symmetry planes.

Thanks,
Eamon

[campolee]

Gerard,

We developed a mesh with the conditions that you recommended above with the two rigid planes of symmetry and applied all necessary boundary conditions (see attached Preview file), but to no avail. The domain again resides in the positive octant, and yet the x-displacement is still negative for some reason. If you could shed some light on this matter, it would be greatly appreciated.

Thanks,
Eamon

[ateshian]

Hi Eamon,

I believe the problem lies with the fact that you constrained all displacement degrees of freedom on the outer surface of the sphere. This makes the problem very stiff, because the initial fluid pressurization produces a nearly isochoric response and the material has nowhere to go upon loading. I checked the paper by Rice and Cleary and they don't impose this boundary condition on the outer surface (only the normal stress and pressure on the inner surface), so I recommend you do that as well.

Another consideration for biphasic problems is that the smallest time step for your analysis should be on the order or greater than h^2/E*k, where h is the thickness of the element where the fluid pressure boundary layer initially develops (i.e., the innermost element in your hollow sphere problem); E=Young's modulus and k=permeability. I haven't checked whether this is satisfied in your input file, but I recommend that you verify it.

I have attached the sample file PressureSliver6.feb I have used to help me answer your question. In my case I use E=1, k=1e-3, and a biased mesh with the thinnest element having h=1e-3. The displacement behaves as expected (you can display the vector plot for displacements in PostView). The nodal 'effective fluid pressure' misbehaves a little in the early time response (t<2) but looks fine after that, up to t=1e4. I think this occurs because of the contact with the symmetry planes. The 'fluid pressure' field, which is evaluated at the element integration points, behaves well throughout the entire time response. You can also see the steady-state response in PressureSliver6SS.feb.

Check it out and see if it helps you fine tune your model until you get satisfactory results.

Best,

Gerard

[siboles]

Hi Folks,

Defining the mat_axis in terms of the local element numbering in this problem is problematic since it has both penta6 and hex8 elements, and the local element numbering will not result in the same coordinate system for all elements. Fortunately, a cylindrical definition for material axis is now possible, so you may want to consider using this rather than making unique materials for the two element types.

http://mrlforums.sci.utah.edu/forums/project.php?issueid=455

Cheers,
Scott