Page 1 of 2 12 LastLast
Results 1 to 10 of 19

Thread: Parameter Optimization for Spherical Indentation of a Viscoelastic Layer

  1. #1
    Join Date
    Dec 2007
    Posts
    906

    Default Parameter Optimization for Spherical Indentation of a Viscoelastic Layer

    The files in this tutorial can be analyzed with FEBio 1.5.2 onward.

    FEBio can be a very effective tool for curve-fitting material properties from experimental measurements. For example, indentation stress-relaxation of a biological tissue with a spherical indenter is a relatively easy experiment to perform, but curve-fitting the resulting load response to an analytical expression is difficult (since analytical solutions are available only for highly specialized conditions, such as Hertzian contact of half-spaces under infinitesimal strains). Since FEBio can solve this contact problem under finite deformation for arbitrary dimensions of the indenter and tissue, it can be used effectively for curve-fitting experimental data.

    To illustrate this process, consider a rigid spherical indenter contacting a cylindrical tissue specimen consisting of a viscoelastic solid whose elastic response is neo-Hookean (SphericalIndenter-axisym.feb). The bottom side of the cylindrical specimen adheres to a rigid substrate. For computational efficiency, since the problem is axisymmetric, only a wedge section is used for the analysis (see the tutorial on Axisymmetric Analyses in FEBio):
    SphericalIndenter-axisym.jpg

    The optimization is performed on three of the four parameters that govern this material: The relaxation time constant t1, the corresponding coefficient g1 in the Prony series of the relaxation function, and Young's modulus E for the elastic part. (Only Poisson's ratio is kept fixed.) The optimization file (SphericalIndenter-axisym-opt.feb) defines these parameters and provides an initial guess, a range (min and max) and a representative scale for each of them. This file also includes the experimental stress-relaxation load response which will be fitted with the optimization analysis. (For this example, the 'experimental' response was generated with FEBio.)

    Note that the experimental load response needs to be adjusted since the actual mesh only uses a thin wedge to represent the true geometry. In this example, the wedge angle is 7.5 degrees, which means that the reaction force predicted from the finite element model is only 1/48 of the actual force (7.5/360). So the experimental reaction force provided in the optimization file needs to be 1/48th of the actual force (in this example), to yield an accurate curve-fitting procedure. Thinner wedges may be used to produce a more faithful axisymmetric analysis, in which case the correction factor for the experimental reaction force should be modified accordingly.

    The command to run this optimization analysis is

    febio -s SphericalIndenter-axisym-opt.feb

    The analysis file SphericalIndenter-axisym.feb must reside in the same directory.

    The optimization analysis for this example converges with 6 major iterations and 28 minor iterations. Each minor iteration takes about 10 s on a Mac laptop, so that the entire optimization analysis completes in just a few minutes. The resulting fit is very faithful,
    SphericalIndenter-axisym-Fz.jpg
    and the fitted parameters (g1=0.9929, t1=99.88, E=1.001) agree very well with the actual parameters used to generate the 'experimental' data (g1=1, t1=100, E=1).

    Gerard

  2. #2
    Join Date
    Jan 2012
    Posts
    2

    Default

    Gerard,

    When will version 1.5.2 be available? Thank you for posting this very useful example!

    Best,

    Kristin

  3. #3
    Join Date
    Nov 2007
    Location
    Park City, Utah
    Posts
    123

    Default

    Hi Kristin - very soon, within a week or two. There are a handful of problems in the Verification Suite that are not producing correct results with 1.5.2. We are going through those problems one-by-one to determine the cause of the errors and correct the code.

    Cheers,

    Jeff
    Jeffrey A. Weiss
    Professor, Department of Biomedical Engineering, University of Utah
    Director, Musculoskeletal Research Laboratories
    jeff.weiss@utah.edu

  4. #4
    Join Date
    Dec 2007
    Posts
    906

    Default

    Hi Kristin,

    FEBio 1.5.2 is now available for download.

    Best,

    Gerard

  5. #5
    Join Date
    Jan 2013
    Posts
    6

    Default

    Hello,

    I am having trouble creating this geometry. From the preview manual example, I can create a ¼ axisymmetric model, but I cannot seem to figure out how to divide the mesh into a smaller wedge. Specifically, when I apply a mesh to the sphere, the only option is to use hexahedrals (which do not form a wedge center like a cylinder does), so I cannot just simply delete away elements to form a wedge. Is there a way to edit the geometry before applying a mesh?

    Thanks!
    Cori

  6. #6
    Join Date
    Dec 2007
    Posts
    906

    Default

    Hi Cori,

    Unfortunately the mesh of the sphere used in this tutorial example cannot be generated in PreView. It was generated in Cubit, exported into Abaqus format and imported into PreView.

    The geometry of the tissue layer was generated with PreView using a Cylinder object meshed using a wedged center. If you don't have (or want) the option of using an alternative meshing software, you can create a hollow sphere in PreView and keep only a single octant. Since the sphere is rigid in this analysis, it won't matter that it spans a larger circumferential angle than the wedge. The computational cost of performing a contact analysis with a sphere octant instead of a wedge will be small.

    Best,

    Gerard

  7. #7
    Join Date
    Jan 2013
    Posts
    6

    Default

    Hello,

    I was using this example to help me create a parameter optimization for spherical indentation of a biphasic layer, and I am coming across a few problems. First of all, when I run the first file, I get an error: “negative jacobian was detected at…” towards the beginning of my simulation. Additionally, towards the end I get a warning: “no force acting on the system.” However, when I look at my results in postview, the stress-relaxation behavior looks similar to what I would expect.

    I also performed the optimization to see what would happen (even with these errors), and it gets through about 4 iterations before it quits and reads: “This application has requested the Runtime to terminate it in an unusual way. Please contact the application’s support team for more information.”

    Do you have any insight as to what is causing these errors?

    Thank you!
    Cori

    SphericalIndenter_CNR_v08.febSphericalIndenter_CNR_v08_opt.feb

  8. #8
    Join Date
    Dec 2007
    Posts
    906

    Default

    Hi Cori,

    First of all, when I run the first file, I get an error: “negative jacobian was detected at…” towards the beginning of my simulation.
    A negative jacobian can occur when the load being applied over that time step is too large, causing some element(s) to invert. FEBio tries to compensate for that by reducing the time step. If the program terminates prematurely it means that the recovery was unsuccessful and you should try to reduce the time step further. If it runs to completion however, the recovery was successful and the solution is valid.

    Additionally, towards the end I get a warning: “no force acting on the system.”
    In biphasic problems, once the response has reached steady state (e.g., in stress-relaxation or creep), there is no further change to the displacement and fluid pressure. FEBio issues this warning because it does not detect any significant change in consecutive time steps. It is a harmless warning message.

    I also performed the optimization to see what would happen (even with these errors), and it gets through about 4 iterations before it quits and reads: “This application has requested the Runtime to terminate it in an unusual way. Please contact the application’s support team for more information.”
    The main reason that a parameter optimization terminates with this type of ominous message is that the analysis for a particular set of material parameters terminated prematurely. This can happen because a negative jacobian was encountered, or the analysis produced a NaN (not-a-number), or if the material parameters assumed values considered unacceptable for that particular material. For example, in your optimization file, you allow Young's modulus to vary from 0 to 5, but 0 is an unacceptable value for a neo-Hookean solid (it will be rejected by that routine's initialization with an error message). This is easily fixed by picking a non-zero lower bound. Similarly, you allow Poisson's ratio to vary from 0 to 0.5, but the exact value of 0.5 is unacceptable (since it produces an infinite bulk modulus, possibly triggering a NaN error), so you should pick a slightly lower value for this upper bound. Keep in mind that in some cases you may pick theoretically valid values but in practice the material might be very soft, possibly leading to negative jacobians. So a certain amount of finagling is needed to formulate a parameter optimization analysis that converges successfully.

    I have a couple of additional remarks based on the files you submitted:

    - You prescribed zero fluid pressure to the entire surface of the tissue, including the region that ends up contacting the indenter. This boundary condition is valid if you assume that the spherical indenter is porous and free-draining. If you intended to have an impermeable indenter instead, you should remove the zero fluid pressure boundary condition from the top surface and employ a biphasic contact interface (sliding2) instead of the facet-to-facet sliding interface you currently use. The biphasic contact interface will automatically enforce zero fluid pressure outside of the contact region and continuity of the fluid flux normal to the contact interface inside the contact.

    - If you do switch to the biphasic contact interface, in case you find poor convergence of the contact analysis, consider switching from BFGS to full-Newton iterations, and from symmetric to non-symmetric storage for the stiffness matrix. These two changes may slow down the speed of the analysis, so use them only if you encounter convergence problems.

    - If you find that the parameter optimization runs a lot of full iterations (e.g., more than 50) while the parameters change by only a tiny amount, consider increasing <f_diff_scale> from 0.001 to 0.01. I have found that this helps in parameter optimizations that involve permeability.

    Let me know if these suggestions improve your results.

    Best,

    Gerard

  9. #9
    Join Date
    Jan 2013
    Posts
    6

    Default

    Hi Dr. Ateshian,

    Thank you for your detailed response (and sorry I am following up so late), it was very helpful! Regarding the negative Jacobian, since the program runs to completion, I did not alter the time steps. Also, I changed my material parameter max and min values so that they are no longer unacceptable for my material. Regarding your additional remarks, I changed the zero pressure condition to a biphasic contact interface. It converges, so I kept it as BFGS iterations and symmetric storage.

    Unfortunately, my optimization is still not working correctly. After approximately 9 iterations, it displayed: “Fatal Error – FEBio error terminated. Parameter optimization cannot continue.” Additionally, the loadcurve I am using to fit the parameters is generated from the original FE model, so I know what values the optimization should output (E=0.5, v=0.3, perm=0.003). However, the optimization parameters were moving farther away from those values before it terminates. Then, I tried narrowing my maximum and minimum parameter values, and the optimization converged on a solution, but it was incorrect. I am unsure how to go about troubleshooting this problem, and I was wondering if you had any further suggestions. Attached are my updated files.

    Thank you!
    Cori

    SphericalIndenter_CNR_v11.feb
    SphericalIndenter_CNR_v11_opt.feb

  10. #10
    Join Date
    Dec 2007
    Posts
    906

    Default

    Hi Cori,

    The problem is with the sign of the reaction force in your optimization file. You need to negate the force values in SphericalIndenter_CNR_v11_opt.feb then you'll get the correct outcome.

    Best,

    Gerard

Posting Permissions

  • You may not post new threads
  • You may not post replies
  • You may not post attachments
  • You may not edit your posts
  •