Modelling cartilage as biphasic-conwise linear elastic material

[lingminl]

Hi,

I'm trying to see if I can use FEBio to reproduce the results from paper "Inhomogeneous cartilage properties enhance superficial interstitial fluid support and frictional properties, but do not provide a homogeneous state of stress." (Krishnan R, Park S, Eckstein F, Ateshian GA. J Biomech Eng. 2003 Oct;125(5):569-77. http://www.ncbi.nlm.nih.gov/pubmed/14618915).

I created a sliver of cylinder to make the cartilage model "axisymmetric", but when comes to the materials section, I'm not sure about it, as I'm kinda of new to biphasic.

My questions are:

1. For just one layer of cartilage (not considering the depth-denpendence yet), how to create the TC nonlinearity according to the given H-A, H+A, lamda2 values in FEBio? What I tried was <solid type="TC nonlinear orthotropic">, where I find it hard to connect the ideas of c1, c2, k, beta, ksi and a, d vectors to the paper given values.

2. how to model the both kr and kz in the permeability type? Should I use <permeability type="perm-ref-trans-iso"> or just constant permeability?

The following is an example of what I tried, the values are not cited from the paper, but only a format as I'm not sure which value I should use:
<material id="3" name="Cartilage" type="biphasic">
<phi0>0</phi0>
<solid type="TC nonlinear orthotropic">
<c1>1</c1>
<c2>0</c2>
<k>100</k>
<beta>4.3,4.3,4.3</beta>
<ksi>4525,4525,4525</ksi>
<mat_axis type="vector">
<a>1,0,0</a>
<d>0,1,0</d>
</mat_axis>
</solid>
<permeability type="perm-ref-trans-iso">
<perm0>0</perm0>
<perm1T>0.00058</perm1T>
<perm1A>0.00186</perm1A>
<perm2T>0</perm2T>
<perm2A>0</perm2A>
<M0>0</M0>
<MT>0</MT>
<MA>0</MA>
<alpha0>0</alpha0>
<alpha1>0</alpha1>
<alpha2>0</alpha2>
</permeability>
</material>

Thanks,

Laura,

[ateshian]

Hi Laura,

1. For just one layer of cartilage (not considering the depth-denpendence yet), how to create the TC nonlinearity according to the given H-A, H+A, lamda2 values in FEBio? What I tried was <solid type="TC nonlinear orthotropic">, where I find it hard to connect the ideas of c1, c2, k, beta, ksi and a, d vectors to the paper given values.

The "TC nonlinear orthotropic" material is an uncoupled material, best suited for nearly incompressible solids. In general, these types of materials should not be used as the solid matrix of a biphasic material, since the porous solid matrix can change in volume as fluid enters or leaves the pore space.

To reproduce the model of the Krishnan paper as closely as possible (keeping in mind that the paper used infinitesimal strain analysis whereas FEBio is a finite deformation code), here is what you need to do for the solid matrix: Use a solid mixture that includes a compressible neo-Hookean solid (type="neo-Hookean") and three orthogonally oriented fibers (type="fiber-exp-pow").

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

You also need to select the orientation of the fibers (parameters "theta" and "phi" in "fiber-exp-pow") to coincide with the radial, circumferential and axial directions of the axisymmetric geometry. The easiest way to do this is to pick a local coordinate system based on node numbering. See the tutorial Holzapfel-Gasser-Ogden material for arteries for a more detailed explanation of how to orient fibers in a local coordinate system. If you created your geometry in PreView using the default coordinate system, and used a "wedged center" meshing method, then add <mat_axis type="local">1,4,5</mat_axis> in your material definition to produces a local coordinate system {x1,x2,x3} aligned with the {circumferential, axial, radial} directions. Then, for each of the three fibers,
circumferential: phi=90, theta=0
axial: phi=90, theta=90
radial: phi=0, theta=0

2. how to model the both kr and kz in the permeability type? Should I use <permeability type="perm-ref-trans-iso"> or just constant permeability?

Because of the above choice of local coordinate system, it is inconvenient to use the transversely isotropic model. Use "perm-ref-ortho" instead, with
perm0 = 0
perm1 = kr, kz, kr
perm2 = 0, 0, 0

So, altogether, your material definition should look like

		<material id="1" name="Cartilage" type="biphasic">
			<mat_axis type="local">1,4,5</mat_axis>
			<phi0>0.2</phi0>
			<solid type="solid mixture">
				<solid type="neo-Hookean">
					<density>1</density>
					<E>0.56</E>
					<v>0.28</v>
				</solid>
				<solid type="fiber-exp-pow">
					<alpha>0</alpha>
					<beta>2</beta>
					<ksi>0.83</ksi>
					<theta>0</theta>
					<phi>90</phi>
				</solid>
				<solid type="fiber-exp-pow">
					<alpha>0</alpha>
					<beta>2</beta>
					<ksi>0.83</ksi>
					<theta>90</theta>
					<phi>90</phi>
				</solid>
				<solid type="fiber-exp-pow">
					<alpha>0</alpha>
					<beta>2</beta>
					<ksi>0.83</ksi>
					<theta>0</theta>
					<phi>0</phi>
				</solid>
			</solid>
			<permeability type="perm-ref-ortho">
				<perm0>0</perm0>
				<perm1>1.9e-3,0.58e-3,1.9e-3</perm1>
				<perm2>0,0,0</perm2>
				<M0>0</M0>
				<alpha0>2</alpha0>
				<M>0,0,0</M>
				<alpha>2,2,2</alpha>
			</permeability>
		</material>

Let me know if this works for you.

Best,
Gerard

[lingminl]

Dr. Ateshian,

Thanks so much for your prompt reply. I applied the material section to the unconfined compression model (see attachment Krishnan2003Axisymmetric.feb) and get it working. The model is a depth-homogeneous model with the femur's layer 1 values, but the loading was only 150N/32 (the cylinder was divided into 32 equal parts) =4.6875N (or actually I should use 150N as the paper mention in the plane strain case here?).

Most of the inputs in your example I understand, except 2 things:

1. In your example, fiber axis is defined right after material line like the following:
<material id="3" name="Cartilage" type="biphasic">
<mat_axis type="local">1,4,5</mat_axis>
<phi0>0.2</phi0>
<solid type="solid mixture">
<solid type="neo-Hookean">...

However, I could not even start running this input in FEBio 1.5.2, so I changed the mat_axis line under the solid type like this:
<material id="3" name="Cartilage" type="biphasic">
<phi0>0.2</phi0>
<solid type="solid mixture">
<mat_axis type="local">1,4,5</mat_axis>
<solid type="neo-Hookean">...

And then it ran smoothly. Was the change appropriate (if so, how can I check the local axis? and was that discrepancy from our different version of input file? I'm using PreView 1.9 to generate the .feb file)

2. In the Krishnan's paper, there is no mention about the phi0 value (or maybe I missed it?). In that case, how and why can we choose this <phi0>0.2</phi0>?

Thanks so much,

Laura,

[lingminl]

Dr. Ateshian,

one more question: when I want to further the model using Inhomogeneous material properties, I chose to simply change their values in each layered elements through a custom-written Matlab code and make these layers still as one object but with 4 different material properties (Krishnan2003Axisymmetric_1.feb).

However, the results are not continuous at the boundary of the layers as I expected, or my expectation is incorrect? Postview picture: postview.jpg PostView plot shear stress: ShearStress.jpg) .

Should I make the the layers individually and create biphasic contact between each layers? Or my current workaround is fine to work but just need some adjustment?

Thanks again,

Laura,

[ateshian]

Hi Laura,

1. The change that you made is correct, it should give you the correct solution.

2. For infinitesimal strain theory with constant permeability, the value of phi0 does not influence the solution. It only become relevant for finite deformation when the permeability is strain-dependent. Using a value of 0.8 is a reasonable value for articular cartilage.

Best,

Gerard

[ateshian]

Hi Laura,

In the Krishnan paper, the material properties measured in the four layers were interpolated using quadratic functions to produce a continuous distribution through the depth. This continuous function was used to evaluate the properties in each layer of elements through the depth. You need to do the same in your analysis in order to reproduce the results from that paper.

Best,

Gerard

[lingminl]

Thanks, Dr. Ateshian. I still have questions regarding reproducing this Krishan 2003 paper's inhomogeneous results:

1. Do I need to change the orientation of the fiber-exp-pow for each depth layer? Or are the depth-depending H+a,H-a, lamda2 all representing the depth-changing fibrils orientations?

2. For the contact between the rigid indenter and the cartilage, should I use the sliding contact or the biphasic contact? If the former, should I make the cartilage articular surface as free draining or not?

3. My model (Krishnan2003AxisymmetricPatella_1.feb) has non-identical effective fluid pressure in circumficial direction (effectivefluidpressure.jpg), was that caused by the not-thin-enough sliver of cylinder?

Thanks again,

Laura,

[ateshian]

Hi Laura,

1. Do I need to change the orientation of the fiber-exp-pow for each depth layer? Or are the depth-depending H+a,H-a, lamda2 all representing the depth-changing fibrils orientations?

No, you don't need to change the orientation with depth. The model used in that earlier paper assumed that fibrils were oriented radially, axially and circumferentially, all directions having the same properties. This was a simplification of the real material behavior of cartilage.

2. For the contact between the rigid indenter and the cartilage, should I use the sliding contact or the biphasic contact? If the former, should I make the cartilage articular surface as free draining or not?

You should use the biphasic contact because it will automatically account for the correct boundary conditions inside and outside the contact region. (Using sliding contact with free draining surface is equivalent to assuming that the indenter is porous and free-draining, which is not the case modeled in that paper.)

My model (Krishnan2003AxisymmetricPatella_1.feb) has non-identical effective fluid pressure in circumficial direction, was that caused by the not-thin-enough sliver of cylinder?

I am not sure about that, but I suspect that you would get better results if the spherical indenter has a finer mesh (same density as the articular surface) and if the indenter is also created as a wedge of the same angle as the articular layer.

I also recommend that you use Newton iterations (instead of BFGS) for biphasic contact problems, and remember to set the matrix to non-symmetric.

Best,

Gerard

[ateshian]

Hi Laura,

1. Do I need to change the orientation of the fiber-exp-pow for each depth layer? Or are the depth-depending H+a,H-a, lamda2 all representing the depth-changing fibrils orientations?

No, you don't need to change the orientation with depth. The model used in that earlier paper assumed that fibrils were oriented radially, axially and circumferentially, all directions having the same properties. This was a simplification of the real material behavior of cartilage.

2. For the contact between the rigid indenter and the cartilage, should I use the sliding contact or the biphasic contact? If the former, should I make the cartilage articular surface as free draining or not?

You should use the biphasic contact because it will automatically account for the correct boundary conditions inside and outside the contact region. (Using sliding contact with free draining surface is equivalent to assuming that the indenter is porous and free-draining, which is not the case modeled in that paper.)

My model (Krishnan2003AxisymmetricPatella_1.feb) has non-identical effective fluid pressure in circumficial direction, was that caused by the not-thin-enough sliver of cylinder?

I am not sure about that, but I suspect that you would get better results if the spherical indenter has a finer mesh (same density as the articular surface) and if the indenter is also created as a wedge of the same angle as the articular layer.

I also recommend that you use Newton iterations (instead of BFGS) for biphasic contact problems, and remember to set the matrix to non-symmetric.

Best,

Gerard

[lingminl]

Thanks so much, Dr. Ateshian. I'll try to modify the model and let you know. Laura,

[lingminl]

Dr. Ateshian,

I've followed the instructions you gave above to re-create the Figure 3a and Figure 4a in Krishnan 2003. So far I got a close plot to them, but still some discrepancies. I've tried to play around the boundary conditions and their corresponding parameters, the contact types and parameters, the initial conditions, the force, the integration parameters, the meshing and densities, the inhomogenous material parameter using different quadratic curve fitting. But still I got 2 main discrepancies that I could not explain:

1. one is the fluid pressure plot: there are larger spikes at the deep zone, which will also create a larger fluid flux in those area, and so even there are fluid going out of the deep zone boundary.

2. one is the 1 principal lagrange strain plot: the there are larger strain values at the superficial surface within the contact.

I'm not sure if those differences were acceptable since "the paper used infinitesimal strain analysis whereas FEBio is a finite deformation code", do you mind take a look at them? Krishnan2003AxisymmetricPatellaThinR100PoroContactFreeDrainingWedgedIndenterNoBiasSymLargeR500In.jpg Krishnan2003AxisymmetricPatellaThinR100PoroContactFreeDrainingWedgedIndenterNoBiasSymLargeR500In.jpg Krishnan2003AxisymmetricPatellaThinR100PoroContactFreeDrainingWedgedIndenterNoBiasSymLargeR500In.jpg

You help or any input is greatly appreciated.

[ateshian]

Hi Laura,

Sorry for the delay in my response but I was traveling extensively.

I think that your results look reasonable in view of the differences between FEBio and the infinitesimal strain finite element code we had written and used in the prior Krishnan study. As you can see from the Krishnan paper, the strains achieved in that analysis approach 0.15 and such values can produce noticeable differences between infinitesimal and finite strain analyses. These differences occur not just because the strain magnitudes cause nonlinear effects, but also because of finite rotations of the fibers. Since a fiber can only sustain tension, and since the fiber strain is evaluated based on its changing orientation, some additional differences may arise with the infinitesimal strain analysis regarding whether fibers are in tension at a particular location through the articular layer.

Best,

Gerard

[lingminl]

Hi Dr. Ateshian,

I really appreciate your time during your busy agenda, and thanks to your reply, which dissolved the final suspicion I had with my modelling.

May I further ask some further questions?

1. In your previous reply

To reproduce the model of the Krishnan paper as closely as possible (keeping in mind that the paper used infinitesimal strain analysis whereas FEBio is a finite deformation code), here is what you need to do for the solid matrix: Use a solid mixture that includes a compressible neo-Hookean solid (type="neo-Hookean") and three orthogonally oriented fibers (type="fiber-exp-pow").

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

I understand how did I get the H-A and lambda2, however, I'm still not sure about how and why alpha, beta and ksi were defined in those values/function. I read through the FEBio user manual 4.1.3.6 "Fiber with exponential-power law" and still not sure. Especially the following comment confused me that why don't we choose beta >2 if we may need transition from compression to tension?

Therefore, use β>2 when a smooth transition in the stress is desired from compression to tension.

2. The conclusion part of your paper "J Biomech Eng. 2007 Jun;129(3):405-12. Equivalence between short-time biphasic and incompressible elastic material responses. Ateshian GA, Ellis BJ, Weiss JA." is very interesting to me, and also I want to apply that to my research:

The equivalence between the instantaneous biphasic and incompressible elastic responses is
valid for any constitutive model, as long as the biphasic constitutive equations reduce to the
incompressible elastic equations when J = 1, as shown for example in Eq. (11).

I want to see if it is possible to create a transversely isotropic transversely homogeneous model without fluid phase that can produce similar results from biphasic TITH model like that in 2003_Krishnan (as I'm going to make the model more complicated in the future, I wan to avoid the biphasic part to save some CPU time). However, since I used compressible solid mixture when reproducting Krishnan's results, my question is, according to 2007_Ateshian, is it appropriate for me to create a solid-phase-only model with incompressible material like Mooney-Rivlin, and compare that with the current biphasic model with compressible TC-NL solid mixture?

Sorry for so many questions, and wish they would not taking you too much time. I understand that my questions may be a little beyond the FEBio alone, and I really appreciate your instruction and patience.

Best regards,

Laura,

[ateshian]

Hi Laura,

I understand how did I get the H-A and lambda2, however, I'm still not sure about how and why alpha, beta and ksi were defined in those values/function. I read through the FEBio user manual 4.1.3.6 "Fiber with exponential-power law" and still not sure. Especially the following comment confused me that why don't we choose beta >2 if we may need transition from compression to tension?

The model used in the Krishnan paper produces a piecewise-linear stress-strain law: Linear in tension and linear in compression, but with different moduli. Therefore the transition from compression to tension is not smooth (the slope changes stepwise). To reproduce this behavior with fiber-exp-pow it is necessary to use beta = 2. alpha is set to zero because it controls the exponential behavior of the fiber stress-strain response, but a linear behavior is needed to reproduce the Krishnan approach. The value for ksi is evaluated based on the theoretical calculation of the fiber elasticity tensor at the strain origin, see Section 5.2.7 in the Theory Manual.

I want to see if it is possible to create a transversely isotropic transversely homogeneous model without fluid phase that can produce similar results from biphasic TITH model like that in 2003_Krishnan (as I'm going to make the model more complicated in the future, I wan to avoid the biphasic part to save some CPU time). However, since I used compressible solid mixture when reproducting Krishnan's results, my question is, according to 2007_Ateshian, is it appropriate for me to create a solid-phase-only model with incompressible material like Mooney-Rivlin, and compare that with the current biphasic model with compressible TC-NL solid mixture?

The short answer to that question is "yes" but it also depends what you want to learn from the finite element analysis. Using an incompressible solid to reproduce the instantaneous biphasic response will produce the same results for the strain (strain_instantaneous_biphasic = strain_incompressible) and the stress (mixture_stress_instantaneous_biphasic = total_stress_incompressible). However, the fluid pressure in the biphasic analysis will not be the same as the Lagrangian multiplier "pressure" of an incompressible elastic solid unless certain equivalence conditions are satisfied as outlined in the 2007 paper.

Best,

Gerard

[yantaaa]

Hi Prof. Ateshian,

I am also thinking about incorporating the CLE T-C nonlinear material to hip joint modelling. As you may have known from my presentation on CMBBE, the model I am using now is isotropic linear elastic on the solid phase, and I would like to move a step further to adopt the CLE T-C nonlinear material. However, I failed to find any paper reporting the hip cartilage properties on its tension modulus and kr.
Do you know any literature on hip tensile properties so that I can give a try in the hip model?
In the study by Cohen et al. 1993 (The influence of transverse isotropy on cartilage indentation behaviour - A study of the human humeral head), the aggregate modulus (or E3 in that paper regarding the near zero of the Poisson's ratio) of the KLM model differs from that of the CLE model. Is this because that the protocol adopted is an indentation model, instead of an unconfined one which I guess may produce a same value of aggregate modulus for CLE and KLM if Poisson's ratio is near zero? Or in other words, would the CLE model curve fitting predictions less depend on the protocol (indentation, confined and unconfined), as compared with the KLM model?

Many thanks.

Regards,
Junyan Li