Modelling cartilage as biphasic-conwise linear elastic material

[ateshian]

Hi Junyan,

For a material such as cartilage I recommend using a continuous fiber distribution instead of the CLE model. See this paper for a justification. For example, for the solid matrix of a biphasic material, you can use a "solid mixture" including a "neo-Hookean" ground matrix and an "ellipsoidal fiber distribution" to capture the tensile response, as described in that paper.

For human hip cartilage a quick PubMed search shows the following two papers that report tensile properties:

Osteoarthritis Cartilage. 2012 Nov;20(11):1268-77. doi: 10.1016/j.joca.2012.07.016. Epub 2012 Jul 31.
Structure-function relationships in osteoarthritic human hip joint articular cartilage.
Mäkelä JT, Huttu MR, Korhonen RK.

Biorheology. 2008;45(3-4):337-44.
Tensile and compressive properties of healthy and osteoarthritic human articular cartilage.
Boschetti F, Peretti GM.
Source

You may have to refit their data to extract the material constants for the model you choose.

Best,

Gerard

[yantaaa]

Hi Prof. Ateshian,

Many thanks for the reply.

As a trial comparision between isotropic linear elastic and Tension-Compression (T-C) nonlinear solid phase for articular cartilage, I adopted the material properties from "Cohen, B., T.R. Gardner, and G.A. Ateshian, The influence of transverse isotropy on cartilage indentation behavior - A study of the human humeral head. Transactions on Orthopaedic Research Society, 1993. 18: p. 185. " for a hip model and a simple indentation model. Results are different between isotropic and T-C, as expected.

Interestingly, for both the simple indentation model and hip model, the fluid support ratio (total fluid pressure on surface / total contact stress*100) for the T-C model is above 100% during early period. I checked the contact stress and fluid pressure on the articulating surface and found that the area of concentrated contact stress has lower fluid pressure than contact stress, but the area away from the peak value region has higher fluid pressure than contact stress. I think this may be because that elements of the eara away from the peak value region are "squeezed" by other elements and thus are in a tension state along the direction vertical to the articulating surface.

Do you think there are other reasons for the >100% fluid support ratio?


So for the T-C model, fluid support ratio may not fully represent the load ratio shared by the fluid / solid. Maybe looking at the stress level for the solid phase would be better.

Is it possible for FEBio/Postview to display solid phase stress?

Regards,
Junyan

[ateshian]

Hi Junyan,

In my experience it is only under special circumstances that the fluid load support magnitude can exceed 100%, typically when the loading regime produces sub-ambient pressure in the cartilage which leads to suction (e.g., under dynamic loading, see this paper).

The reason you are getting (slightly) greater than 100% could be for a variety of reasons which I cannot figure out without taking a look at the model and understanding exactly how you are calculating the fluid load support for that problem. The more information you can provide me the better I can help answer your question.

Currently PostView does not provide the effective (solid) stress sige as a direct output. You can export the stress sig and fluid pressure p to a file and you can calculate the effective stress as sige=sig+p (only for the normal components of the stress, since the total and effective shear components are the same). You can load this data back into PostView to display it (click on the Data tab and select Add->From File...).

Best,

Gerard

[yantaaa]

Hi Prof. Ateshian,

Thanks again for the reply.

Attached is the simplified model replicating the situation of a hip joint.
Please have a look.

For the results, "squeezed" elements can also be detected. For example, element 79, 82 and 85 are squeezed along the tangent direction to the articulating surface and are in a tension state along the direction vertical to articulating surface. So for these elements, the fluid pressure should be higher than contact stress, and the fluid support ratio (if calculated as fluid pressure/contact stress) for these element should be above 100%.

I calculated the fluid support ratio by selecting the articulating surface and "integrate" the fluid pressure, contact pressure/traction and also 3rd stress. However, I think 3rd stress may not represent contact stress very well for all the elements on the whole articulating surface.
Is there a more accurate way to calculate fluid support ratio to consider the elements of >100% fluid support ratio as =100% ratio?

For loading data back to Postview, I added this line <element_data data="p" name="p" file="Sp2.dat"></element_data>. But the fluid pressure loaded by this file cannot be displayed in Postview. Did I conduct this not correctly?

Many thanks.

Regards,
Junyan

[ateshian]

Hi Junyan,

The model you attached uses a very coarse mesh and I am not confident that the results from that analysis are meaningful. I regenerated a similar geometry but used a finer mesh with 10 elements through the layer thickness and mesh biasing to produce a finer mesh near the contact surface and the interface with the rigid substrate, please see the attached files simple_spherical3.feb and simple_spherical4.feb.

To help evaluate the fluid load support I added <var type="contact force"/> and <var type="fluid force"/> to the output plotfile variables, which provide the net contact force components and the net fluid force components on the contact surfaces. The contact force is the integral of the contact traction over the contact surface; on each face the contact force is directed normal to that face since the contact interface is frictionless. Similarly for the fluid force, which is evaluated using only the fluid pressure contribution to the contact traction. Since the displacement is applied along the y-direction, the fluid load support may be calculated from the ratio of the Y-component of the fluid force to the Y-component of the total force.

In the file simple_spherical4.feb the solid matrix only contains a neo-Hookean solid. For this file the fluid load support at the last time point is 37.1/37.55 = 98.8%, which is less than 100% as expected. In the file simple_spherical3.feb the solid matrix also includes a spherical fiber distribution (a continuous fiber distribution which is isotropic). Because the fibers increase the tensile stiffness considerably, the contact load produced in this highly congruent contact configuration rises to much higher values (and the solution terminates a little prematurely). The fluid load support at the last time point is 1199/1191 = 100.6 %. I think that this result simply reflects numerical round-off error (since the fluid pressure is evaluated at the nodes but the contact traction is evaluated at integration points) and that one should interpret the outcome to mean that the fluid load support is effectively 100%. Indeed, in our earlier theoretical analysis (see eq. 27) we found that the fluid load support would approach 100% as the ratio of tensile to compressive stiffness increased to infinity. I suspect that refining the mesh further (along directions normal and tangential to the contact surface) would reduce this round-off error. Further evidence in support of this interpretation is the observation that extracting the Y-fluid force and Y-contact force from the opposite contact surface produces 1199/1197 = 100.2%.

Please try these out and verify those results for yourself.

Best,

Gerard

p.s. To import data back into PostView you must put it in the correct format, see this post.

[yantaaa]

Hi Prof. Ateshian,

Thanks a lot for the detailed reply.

I agree that the round-off error may be a reason for the above 100% fluid support ratio. However, I still think that the >100% value is much related with the "squeezed elements" effect.
Around the contact region (near the axis of symmetry in simple_spherical3.feb) of the articulating surfaces, elements of one layer of cartilage are compressed by the opposite cartilage. Under such circumstances, the solid phase helps the fluid phase to resist the compressive force vertical to the articulating surface, and the fluid pressure is lower than the contact pressure/traction for this region.
However, for the elements near the edge of contact area or outside the contact area, the compressive effect caused by the opposite cartilage layer may not be as high as the squeeze effect by the nearby elements nearer to the contact centre (axis of symmetry). As a result, these elements may be in a tension state vertical to the articulating surface and the solid phase helps resist such tensile stress. In this case, the fluid pressure is higher than the contact pressure.

Please have a look at the figure below. I integrated the fluid pressure/effective fluid pressure and the contact pressure/traction for two regions of the articulating surface respectively: the half region near the axis of symmetry and the other half. For the "other half" region, as highlighted in the figure, the integrated fluid pressure is 127 and the integrated contact pressure is 111. So the fluid support ratio is >110% for this region, which is also the case for effective fluid pressure or contact traction.

If this is true, maybe it is necessary to calculate the fluid support ratio in another way, because fluid support ratio for such problems should not be above 100%. For elements with higher fluid pressure, we need to make its fluid support ratio as 100%, as the higher fluid pressure is not for load supporting any more.

Regards,
Junyan


fig.jpg

[ateshian]

Hi Junyan,

For a given element (or sub-region) in the contact interface it is indeed possible for the fluid pressure magnitude to exceed the magnitude of the total contact stress. This can happen as a result of the local condition for a particular element, such as the squeezing effect you mention. The definition of fluid load support I have been using is evaluated from the integral of the fluid pressure over the contact area, divided by the integral of the contact pressure over the contact area. It is this net effect that I don't normally expect to exceed 100%, but keep in mind that this is not a general theoretical requirement as mentioned in my first response to your question. So it is ok to find that the fluid pressure magnitude exceeds the contact stress magnitude in some locations on the contact interface. The total normal contact traction is t = -p + te, where p = fluid pressure and te=effective normal contact traction (resulting from the solid matrix strain). Mathematically (and physically) it is possible to have |t| < |p|, |t|=|p| or |t|>|p|. For a contact interface however, keep in mind that t ≤ 0. When p > 0, the ratio p/|t| exceeds 100% only when te>0 (tensile effective traction).

Best,
Gerard

[Squid]

Dear All,

I am using a solid mixture that includes a compressible neo-Hookean solid (type="neo-Hookean") and three orthogonally oriented fibers (type="fiber-exp-pow") to model the top layer of a two layer indentation model. The model converge but something really strange happens. Sometimes FEBio crashes at the starting of the analysis. However, strangely, without applying any modification to the .feb, after one or two trials trying to launch the simulation, FEBio runs regularly. How could I do to avoid this problem? It does not happen using a different constitutive law (just changing the keywords related to the material model typology).

Thanks a lot

Regards,

Leo

[ateshian]

Hi Leo,

In FEBio, fiber models such as "fiber-exp-pow" can only sustain tension. Their modulus only enters into the stiffness matrix calculation when the fiber strain is positive. This means that at the strain origin (at the start of any analysis), there is a discontinuity in the stiffness matrix that causes problems with convergence. This is an inherent behavior of tension-only materials. It may be easier to overcome the convergence problems at the start of an analysis by ramping up the prescribed displacement or load (instead of prescribing a step displacement or load) and using a smaller initial time increment.

Best,

Gerard

[Squid]

Dear Prof. Ateshian,

thanks a lot for the answer. In this case the anomaly that I am facing is generated really at the beginning, when I launch FEbio. The application crashes before the simulation starts. It is strange that for the same file the crash happens randomly. I have attached the file to give you a better idea. When I launch the file there is the probability that the application crashes: "FEBio has stopped working". If you run the file that is attached, if it crashes, wait and retry. Generally after two or three attempts FEBio runs and simulation does converge.
Please let me know if you need further details. What it is strange is that the simulation converge...the problem seems arising before the attempt to generate the stiffness matrix.

Thanks a lot

Regards,

Leo

[ateshian]

Hi Leo,

The problem appears to be with the <mat_axis type="local"> statement in that file. You have provided only two node numbers, but this field expects three local node numbers (or none). FEBio does not check if the user has entered three values. As a result, the third value is automatically set to 0 and this appears to cause an error (but not consistently).

Also please note that the node numbers that should be provided in this field should vary between 1 and the number of nodes in that element (e.g., 1 to 8 for a hexahedral element) since you have selected local node numbering. But your file uses 1718 and 1719, which seem to refer to absolute node numbers. These choices are incorrect and need to be fixed.

Best,

Gerard

[Squid]

Dear Prof. Ateshian,

thanks a lot for the answer. I see, and how I could define the coordinate system to align the fibers tangentially to the radial direction? I created the mesh in Abaqus, how should I set the node values? Can I always use <mat_axis type="local">1,4,5</mat_axis> ?

Thanks a lot.

Regards,

Leo

[Squid]

Dear Prof. Ateshian,

rebuilding the mesh with Preview the fibers direction are assigned correctly. In this way I can use:
<mat_axis type="local">0,0,0</mat_axis>
...
<theta>0</theta>
<phi>0</phi>

to orient them in the radial direction. Thanks a lot for the hints.

Regards,

Leo

[ateshian]

Hi Leo,

I am glad you could resolve that issue.

Best,

Gerard

[Squid]

Dear Prof. Ateshian,

considering a continuous fiber distribution and a "Holmes-Mow" solid, are these relations always valid?:

The properties of a "neo-Hookean" solid are "E" (Young's modulus) and "v" (Poisson's ratio). They are related to H-A and lambda2 as follows:
H-A = (1-v)E/(1-2v)(1+v)
lambda2 = v E/(1-2v)(1+v)

The properties of all three "fiber-exp-pow" solids are "ksi", "alpha" and "beta". To reproduce the Krishnan paper, you need to select
alpha = 0
beta = 2
ksi = (H+A - H-A)/4

and where can I find references to justify them?

Thanks a lot

Regards,

Leo