Welcome to the FEBio Forum. This forum serves the community of FEBio and FEBio Studio users and developers. Forum participants are encouraged to post questions, as well as answer posts from others. A broad level of participation is encouraged, to create a vibrant community that helps improve the quality and usefulness of these open-source/free software products. Moderators are here to assist with explaining novel features, addressing bug reports and reviewing feature requests, but the effectiveness of the forum depends critically on the participation of experienced users who can assist novices or share ideas and models that explore challenging problems. Please feel free to join in this effort. You can subscribe to forums by pressing the "Subscribe" button at the top of the forum. This will allow you to stay up to date on recent activity on the forum.
The FEBio software downloads and knowledgebase can be found here.
Need help in modeling indentation stress relaxation of cartilage
I did try playing around with the scale factor and the upper and lower bounds but now the model has problems converging. I have attached the optimization file for you to review. Please let me know what I can do to overcome the problem.
I was wondering if you got a chance to look over the files that I last sent you. I have been trying to vary the scale factors, min, max values and having no luck in finding a solution. Could you please suggest something that I should try?
Since the analysis of this model was taking too long to run, I was not able to dedicate time to figure out what might be the problem. To speed it up, I re-created the model with a coarser mesh (SFDCartInd6mNcoarse.feb). Then, to help me figure out what might be good initial guesses for the optimization, I ran a few cases and adjusted the properties manually to get a reasonably similar response to the experiment (in particular, I tried to match the peak load at the end of the ramp and the steady-state load). That required a half-dozen runs and I found that E=20, ksi=3 and k=2e-5 might be reasonable initial guesses. Then I started the optimization with those guesses and I constrained the range of k to be within factors of 0.1 and 10 of this initial guess. For E and ksi I estimated ranges using factors of approximately 0.25 to 4 of the initial guesses (OptD02LCLatestmNcoarse.feb). These are just suggestions of how to approach an optimization problem when the material parameters are not known. At the end of the optimization (which took a few minutes to run) the fit looked good as shown in the attached figure, producing R^2=0.86 (coefficient of determination). The final fitted parameters were E=30.78, ksi=1.833 and k=5.330e-5. SFDCartInd6mNcoarse.png
For completeness, you now need to check mesh convergence (i.e., that the curve-fitting produces substantially similar values when using a finer mesh, e.g., within 5% or 10% of the coarser mesh). If the results are significantly different, the mesh I used for the above example is too coarse and you would need to make it finer.
First of all thanks a lot for taking out time from your busy schedule to work on the model. I greatly appreciate it.
The model converges on a finer mesh as well with the following final values
E = 30.206
k = 4.55e-05 and
ksi = 1.724.
The fit looks good (please see attached figure) with RSQ value of 0.96. Do you think the model can be further improved? I am thinking about using it for completing the analysis on rest of the samples unless you say otherwise.
If you would be kind enough to send me your prv or even explain to me how you made the symmetry plane (part 3) for the initial geometries of the model?
We have recently introduced a "symmetry plane" option under Physics->Add contact... in PreView. You should be able to select the rotated faces of a wedge geometry to this symmetry plane and prescribe a penalty parameter (a large value, e.g. 1e3 or 1e5). Please try it out and let me know if that doesn't work for you.
Thats me added a symmetry plane. When I run it I am getting negative Jacobian detected. Any idea as to why? if you have time to look iv attached my .feb file bellow. Many thanks.
You are getting negative jacobians for the following reasons:
1) You are running a load-control analysis but there is no initial overlap between the indenter and the biphasic layer. So, at the very first iteration of the first time point, there is no force from the cartilage layer countering the one applied to the indenter. The indenter moves to infinity, which produces negative jacobians. In your model you need to translate the indenter down so that it overlaps a little bit with the cartilage layer (e.g., translate along Z by -0.002).
2) The load you applied on the indenter is positive, which means the indenter will move up instead of down. You need to negate the load value (or the loadcurve).
3) The load magnitude is probably too large for the properties you are using for this problem. Remember that a wedge only occupies a portion of the circumference (e.g., a 3 degree wedge represents 120th of the full circular geometry), so the load you must apply should be scaled down accordingly (e.g., by a factor of 120).
Ok great, that all makes sense. Iv now got the tutorial working. Thanks.
Another question. What would be the most effective way to change parameters such as Youngs Modulus, permeability and ksi and to study these changes to the overall model? What changes happen over time and what's different?
Im not sure if that makes sense but I hope it does.
Young's modulus controls the equilibrium response to a creep problem (such as the example you are analyzing): It determines the relation between the applied (constant) load and the equilibrium creep deformation. So, under a constant load, increasing E decreases the creep deformation.
The permeability and ksi control the transient creep response (how long it takes to reach equilibrium, and the shape of the transient response). More specifically, the characteristic time constant for reaching equilibrium is proportional to the square of a characteristic measure of length (e.g., the indenter radius or the cartilage thickness) divided by the product of ksi and permeability.
Magic, that helps explain what I was observing in post view.
If you have time could you please take a look at my model (iv attached) and tell me what i'm doing wrong in terms of loads. Iv made sure to negate the load curve of the indenter and also to have an initial overlap of -1e-03 in the z axis. But it still won't show any deformation when running the analysis. I also found it difficult to get the sphere to form a wedge centre but does this even matter if its rigid?
The load you applied does not cause too much deformation. I was able to increase it by a factor of 10 and also by a factor of 100 and in both cases the analyses converged and the solutions shows more significant deformation.
When you choose a non-symmetric stiffness matrix, you should switch the quasi-Newton method from BFGS to Broyden. But the best results for contact analyses is achieved using full-Newton iterations (set max_ups=0).
I question about symmetry. I have made two spherical indentation models. One model is 1/4 tissue layer and 1/4 indenter. The other model is a 3 degree tissue layer and 1/4 indenter. The idea behind using the 3 degree wedge was reduced computational time. Will my output analysis results be identical in both models? Ideally I would like to make the spherical indenter also 3 degrees but I have trouble doing this with the editable mesh. Or does it not matter how much of the indenter is present as its a solid?
In principle (i.e., if the penalty parameter for the symmetry plane is high enough), the axisymmetric analysis results between a 1/4 (90-degree) model and a narrow wedge (3-degree) model should be the same for all the nodal and element variables (pressure, displacement, stress, etc.). The indenter load response will differ by a factor of (90 degrees / 3 degrees) = 30.
It is okay if the mesh of the rigid spherical indenter spills over the 3-degree wedge, that should not affect the analysis results.
Tweet
Comment