Dynamic analysis problems

[YIDAN]

Hi,

I'm using a knee model to simulate the 10NM valgus torque with prestrain. I defined one FEBio step from 0 to 3 sec, including the prestrain application step(from 0 to 1s), flexion step(from 1 to 2s), and loading step(from 2 to 3s). For boundary conditions, the femur and tibia are rigid bodies, the tibia is fixed, and the femur is free except the flexion rotation. The analysis type in control section is "static," and it's the "full-Newton" method. When the flexion angle is 0deg/30deg/60deg, the simulations work well, and results look good to me. However, when the flexion angle is 90 deg, the simulation keeps running for more than 2 days (the computer has 4 cores/8 threads, 16GB RAM) and only finished 90% of the flexion step after 2 days(still running, no error termination). It had more than 600 iterations in each time step. So I tried to change the analysis type from "static" to "dynamic." Then, the simulation results look weird after using the dynamic type: the fiber stretch of the end of the ligaments attached to the femur is super large, which is approximate 5 times than usual. Two ends of the ligaments use "rigid nodes" attached to the femur and tibia, respectively. I was wondering why this happened? How can I solve it? Should I use other analysis types instead of the dynamic? Or should I use the Tied contact to define the lig-bone-lig?

Any suggestions/help, I will appreciate it a lot!!

Best,
Yidan

[ateshian]

Hi Yidan,

I assume that your knee model involves contact interfaces. FEBio's sliding contact interfaces are not yet designed to work with dynamic analyses. That's probably why you are seeing unreasonable results.

In general contact interfaces only satisfy the momentum balance equations. To run properly under dynamic analyses they also need to satisfy the energy balance. We plan to implement that at some point in the future, but there is no specific time frame for completion of this task as of now.

Best,

Gerard

[apaz]

Hi Professor Ateshian,

Based on your answer, I wonder some questions. Doesn't FEBio have any contact definition for dynamic problems? Can Yidan use FEBio to solve the problem is trying to model? Is it possible to guess how much would differ answers between static and dynamic? How does it affect when the materials in contact are biphasic or viscoelastic?

Regards,
Alex

[ateshian]

Hi Alex,

Doesn't FEBio have any contact definition for dynamic problems?

None that conserve energy, currently.

Is it possible to guess how much would differ answers between static and dynamic?

The outcome of dynamic analyses using contact interfaces may be unpredictable. For example, if a cube of elastic material is dropped under the action of gravity to make contact with a rigid plane, the cube should normally bounce back to the initial height it was released, over any number of cycles. However, running that analysis in FEBio shows that for some cycles the cube bounces back to a lower height while for other cycles it bounces back much higher (energy added or lost with no obvious consistency). On the other hand, there may be cases when dynamic contact analyses give seemingly reasonable results (there are some examples of those problems in the test suite). For example, if energy is consistently lost but not gained, the results may be acceptable (as though the system has damping, but it would be numerical damping in this case). You can try to tweak the time integration (Newmark) parameters used in dynamic analyses to see if you can promote consistent numerical damping. The Newmark parameters are denoted by <beta> and <gamma> and they can be specified in the <Control> section by also letting <rhoi>-2</rhoi>. See section 6.1 of the Theory Manual for some more details about Newmark integration. There is also the parameter <alpha> which is needed to provide a complete time integration scheme. See section 6.2.3 of the Theory Manual for a description of the Generalize alpha-method to understand the meaning of alpha. (Using the spectral radius rhoi in the range 0 to 1 automatically employs the Generalize alpha-method, which sets the values of alpha, beta and gamma given a single value for rhoi, so keep that in mind when you choose to tweak alpha, beta and gamma on your own.)

How does it affect when the materials in contact are biphasic or viscoelastic?

The biphasic materials in FEBio cannot accommodate dynamic analyses. To do dynamics with biphasic materials it is necessary to add the fluid velocity as nodal degrees of freedom. The biphasic material currently uses only solid displacement and fluid pressure as nodal degrees of freedom. We are currently developing a biphasic-FSI formulation that accommodates dynamics, but it won't be available for another few months. Even in that case, we would not release contact interfaces for this type of domain at the same time, that would require a lot of extra work (so more time).

In FEBio3, energy-conserving dynamic analyses are well-developed for elastic and rigid materials (but no contact). To conserve energy in elastic materials, we need to use a special scheme that requires the evaluation of the elastic material's strain energy density at the current and previous time steps. To enforce energy conservation in these analyses, it is necessary to set the spectral radius rhoi to 1. For standard viscoelastic materials, there is no method to evaluate their strain energy density (and free energy cannot be conserved). So using rhoi=1 will not necessarily improve the analysis. Nevertheless, viscoelastic materials can be used in dynamic analyses (along with elastic and rigid materials).

Best,

Gerard

[apaz]

Hi Professor Ateshian,

Thank you for this elaborated response, I appreciate you share your knowledge and time. It is starting to be clear to me. But, why in works like those of Li et. al. (https://journals.sagepub.com/doi/10....54411914537617), the models work well? Maybe I am misunderstanding some definitions.

Regards,
Alex

[ateshian]

Hi Alex,

That study (and others that use biphasic materials in FEBio) use quasi-static analyses (which is the default analysis type for biphasic materials), so they can model transient responses but don't account for inertia terms in the solid and fluid constituents. In practice, unless you are interested in wave propagation analyses in a deformable continuum, you don't need to use a dynamic analysis. (It is different for rigid bodies, where many common analyses may require dynamics, or cardiovascular fluid mechanics where inertia terms are not negligible, especially in the fluid). For elastic and biphasic materials the rule of thumb is as follows: If the characteristic time of loading is much larger (e.g., at least ten times larger) than the time it takes for a wave to propagate through the deformable medium, you don't need to do a dynamic analysis. In biological tissues which are hydrated, the speed of sound in water is approximately 1480 m/s. If you are dealing with cartilage layers whose thickness is on the order of ~ 1 mm, it takes less than 1 microsecond for a dilatational wave to propagate through the tissue. As long as the characteristic time of loading is greater than ~ 10 µs, it is not necessary to perform a dynamic analysis for cartilage.

This is not the same as saying that dynamics is not important for locomotion studies. One can determine the joint forces using a dynamic analysis of locomotion, but when those forces are applied to your finite element model of the knee, you don't need to include inertia terms.

Best,

Gerard

[apaz]

Hi Professor Ateshian,

Thank you, one more time, for the clarification. Now, it is completely clear to me.

Regards,
Alex