Axisymmetric analysis of biphasic indentation
This tutorial examines the indentation of a biphasic layer with a rigid impermeable frictionless indenter under creep loading. The purpose of this tutorial is to demonstrate the 3D modeling of axisymmetric problems using wedge geometries and a plane of symmetry not aligned with the coordinate planes. The tutorial also shows how to create biased meshes, and attaching and welding congruent meshes.
This is an intermediate tutorial and assumes some familiarity with creating FEBio models in FEBio Studio. If you are new to FEBio and FEBio Studio, you may want to follow the Intro to FEBio Studio tutorials first.
To start the tutorial, open FEBio Studio and start a new model (menu File\New Model). Select the Biphasic Analysis option in the New Model dialog box and click OK. Then, follow the steps below.
Step 1: Creating the geometry
Since FEBio does not provide 2D axisymmetric elements, an axisymmetric analysis may be simulated in 3D using only a wedge (slice) out of the true 3D geometry. Reducing the wedge angle increases the faithfulness of the 3D analysis to the 2D axisymmetric analysis. In this tutorial we will use a wedge angle of 3 degrees. The first objective is to create the wedge representing the biphasic layer and mesh it as shown in the figure below.
Note that the mesh is biased along the radial and axial directions. The axial direction has a dual bias (the mesh is finer at both ends and coarser in the middle). This type of biasing is a standard feature that can be implemented by simply checking a box and entering a bias value as shown below. The radial direction also shows a dual bias, but the location of the dividing line is not at the midpoint of the radial extent. To reproduce this non-standard bias, two geometries must be created, meshed congruently, then attached and welded.
First, create the inner Slice with Radius = 1, Height = 1, Angle (deg) = 3. Second, create a Tube with Inner Radius = 1, Outer Radius = 3, Height = 1. Mesh the slice using Slices= 1, Segments = 12, Stacks = 16, Z-bias = 1.4, R-bias = 0.8, Z-mirrored bias = checked. Mesh the tube to be congruent with the slice along the axial and circumferential directions: Slices = 30 (to get 3 degrees per slice), Segments = 12, Stacks = 16, z-bias = 1.4, r-bias = 1.4, z-mirrored bias = checked. The mesh should now look as shown in the figure below.
Convert both geometries to editable meshes. For the outer tube, delete all elements except the outer slice that meets the inner slice (next figure).
Next, we’ll merge the two objects into a single object. To do this, select both objects. Then, select the menu Edit ∖ Merge objects. In the dialog box that appears make sure the weld option is checked and press OK. This completes the geometry and mesh for the tissue layer.
To create the indenter, activate the Create panel again and create another slice with Radius= 1, Height= 1, Angle (deg)= 3, at the position x= 0, y= 0, z= 1. Mesh it with Slices= 1, Segments= 1, Stacks= 1.
For a creep analysis, where a load is prescribed on the indenter, the contact between the indenter and tissue layer will fail unless the two contacting surfaces overlap at the initial time point. Therefore, a small initial overlap should be created between the indenter and the tissue layer. Make sure the indenter is selected and select the menu Edit ∖ Transform. In the dialog box enter the value -1e-4 in the Relative Z field and press OK.
Step 2: Setting up the materials
Create a Biphasic material (call it Tissue) with a neo-Hookean solid and constant permeability. Assign this material to the tissue layer. Set the properties as follows:
- solid volume fraction = 0.2
- Young’s modulus = 0.4
- Poisson’s ratio = 0.02
- permeability = 0.02
Create a Rigid Body (call it Indenter) and assign it to the indenter.
Step 3: Setting up the boundary conditions
To enforce an axisymmetric response, a symmetry plane should be created which encompasses the deformable regions. Select the surfaces of the biphasic material which lie on the rotated wedge faces. Then select Physics/Add Surface Constraint and select the Symmetry Plane option. Set its parameters as follows:
- laugon = off
- penalty factor = 1e6
Fix the displacement of its front faces along y. Also fix the x and y displacements of the edge coinciding with the axis symmetry. Fix the x, y and z components of the bottom face of the tissue layer.
Add a Rigid Constraint to the Indenter, with Fixed Displacement/Rotation enforced for x, y, Rx, Ry, Rz, leaving z free.
Add a Rigid Constraint to the Indenter, with Prescribed rigid force enforced for z and value of 1.0. Set the load curve for this force to Curve Type == Step with a value of -0.002 at t== 1. Note that the actual force is 120 times this value, since we are using a 3-degree wedge (360/3== 120) in this axisymmetric analysis.
Add a sliding-biphasic contact interface between the top surface (inner and outer slices) of the Tissue (Primary surface) and the bottom surface of the Indenter (Secondary surface). You may need to hide the Indenter to select the Tissue contact surface underneath it, and vice-versa to select the Indenter contact surface. Use the following non-default parameter:
- Enforcement method = AUGLAG (augmented Lagrangian)
- auto-penalty = Yes
- penalty = 1
- symmetric stiffness = off
Step 4: Defining the analysis step
Add a Biphasic step and use the following non-default parameters:
- Analysis: TRANSIENT
- Time Stepping
- Time steps = 10000
- Auto-time stepper
- optimal iterations = 50
- Nonlinear Solver
- Max updates = 0
- Linear Solver
- Matrix storage = Non-symmetric
Add a load curve to the dtmax time stepper parameter. Edit the loadcurve to be of Curve Type= Linear and include the points (0,0) and (1000,200).
Run the analysis. The resulting indenter creep displacement curve and fluid pressure distribution are shown below.
