Biphasic simulation inconsistencies with Ogden material

[lefeverj]

I am continuing to work on a biphasic, pressurized cavity problem with the Ogden material. The spherical geometry consists of two hollow spherical subregions with biphasic contacts to enable local refinement in mesh density. The simulations I have been working on lately have terminated early due to a negative Jacobian.

Repeated runs of this simulation on the PC (with exactly the same input file) consistently yielded the same failure in the same element and time step. However, repeated runs of this simulation on OS X (again, with exactly the same input file) yielded a negative Jacobian in different elements at different time steps (and different from the element and time step on the PC). All failures occurred in elements along the x-axis on the inner subregion, though the particular element varied from run to run. Attached are the feb file and three log files from various runs of the simulation (one on the PC and two on OS X).

In practice, I have done this simulation on a model similar to the one attached but with three subregions and a higher mesh density. In those simulations, the failures occurred more randomly, in both the innermost and outermost subregions and not always along edges or surfaces. The failures in the outer subregion typically occurred in elements that had not deformed significantly from their unloaded state.

As an aside, the problem usually does not occur at all when the Ogden material is replaced with an isotropic elastic material definition.

Can anyone help to solve this problem?

Thanks,
Joel

[ateshian]

Joel,

Ogden's constitutive relation was formulated for incompressible elastic solids. In FEBio it is implemented as an uncoupled material to help enforce incompressibility. In principle, it should always be used with a bulk modulus k which is much larger than the coefficients c[i]. Setting the bulk modulus to a low value (e.g., a value comparable to c[i]) invalidates the analysis because the principal stretches used for evaluating the stress tensor are formulated under the assumption that the deformation is isochoric.

In contrast, the solid matrix of biphasic materials is expected to be compressible because the pores of a porous-permeable matrix must be able to change in volume as fluid enters or leaves the matrix. This means that using an Ogden material for the solid matrix is an inconsistent combination of two frameworks, because Ogden is incompressible while the porous solid matrix must be compressible. For these reasons, FEBio's implementation of biphasic materials does not account for the possibility that the solid matrix is (nearly-) incompressible. I believe that's the reason you're not getting good convergence when you use the Ogden model, whereas using a compressible elastic matrix produces acceptable results.

You can double check this issue by analyzing two alternatives:

(1) Perform your analysis for a non-biphasic (purely elastic) problem using the Ogden material (with an elevated bulk modulus). In principle this problem should run well.

(2) Perform your analysis with a biphasic material whose solid matrix is a Mooney-Rivlin solid (with an elevated bulk modulus). I expect that this problem will not run well, for the same reason that Ogden fails.

If you do try these alternatives, please share your results.

Is there a compelling reason that you should use an Ogden model for the solid matrix of a biphasic material?

Best regards,

Gerard

[smithjh]

Gerard,

My research group, of which Joel is a member, is interested in modeling brain tissue as a biphasic material. Experimental characterization of brain tissue is rather limited, with the two primary works being

-- Miller and Chinzei, Mechanical properties of brain tissue in tension, Journal of Biomechanics, 35:483--490, 2002

-- Franceschini et al., Brain tissue deforms similarly to filled elastomers and follows consolidation theory, Journal of the Mechanics and Physics of Solids, 54:2592--2620, 2006.

In both of these studies, the behavior of brain tissue under large deformation is characterized using the incompressible Ogden material model. Therefore, it is our desire to model brain tissue as a biphasic material with a (compressible or incompressible) Ogden material model for the solid phase.

I would like to follow-up on your statement that "the solid matrix of biphasic materials is expected to be compressible because the pores of a porous-permeable matrix must be able to change in volume as fluid enters or leaves the matrix." Is this a general rule for biphasic material or just a requirement because of how FEBio is implemented for biphasic materials? Is it not possible, for example, that the bulk volume of a biphasic material increases through fluid absorption and expansion of the pores but that the solid phase be nearly incompressible and not change volume during the process? Could the cells not lengthen and thin to allow for larger pores while themselves maintaining the same volume?

If it is not possible (in general or just for FEBio) to have an incompressible solid model within a biphasic material, would it be possible for a compressible version of the Ogden material model be added to FEBio?

As an aside, I must admit that we are confused by the consistency of your response with a forum posting by Jeff in which he said that "'incompressible' is a mis-label" for the Ogden material and that it "should be fine for compressible values of the bulk modulus."

Finally, Joel re-ran your suggested alternatives. His results were:

1. For a monophasic simulation with the Ogden material using the original bulk modulus or a value 100x larger, this converged, as you expected. Additionally, for a biphasic simulation with the Ogden material using a bulk modulus 100x larger, this still failed to converged, as I think you would expect.

2. For a biphasic simulation with the Mooney-Rivlin material using the original bulk modulus or a value 100x larger, this actually converged, in contrast to your expectations.

Any further insight into this issue would be greatly appreciated.

Thanks,
Josh

[smithjh]

Gerard,

In reviewing the FEBio Theory Manual, I have determined that my misunderstanding came from confusing "solid phase" or "solid constituent" with "solid matrix." I now understand that the solid constituent is incompressible and the solid matrix must be compressible. I further understand why the incompressible Ogden material model is inconsistent with a biphasic formulation.

However, I would still like to know if you think it is possible (with reasonable effort) for a compressible Ogden material model to be included in FEBio. Considering that the aforementioned experimental studies for brain tissue, this would be greatly appreciated.

Thanks,
Josh

[weiss]

Joel,

Ogden's constitutive relation was formulated for incompressible elastic solids. In FEBio it is implemented as an uncoupled material to help enforce incompressibility. In principle, it should always be used with a bulk modulus k which is much larger than the coefficients c[i]. Setting the bulk modulus to a low value (e.g., a value comparable to c[i]) invalidates the analysis because the principal stretches used for evaluating the stress tensor are formulated under the assumption that the deformation is isochoric.

Hi Gerard - I just wanted to comment that the Ogden material should work fine with a relatively low bulk modulus for an elastic analysis. The constitutive model is uncoupled and uses the three-field "Fbar" element, so the principal stretches in the model are the deviatoric principal stretches.

The constitutive model is described in Simo and Taylor's paper, "Quasi-incompressible finite elasticity in principal stretches: Continuum basis and numerical algorithms. Computer Methods in Applied Mechanics and Engineering, 85:153-173, 1991.

The model can be used with the augmented Lagrangian method to enforce incompressibility to a user-defined tolerance, but the user must still supply an "initial" bulk modulus and then turn on the augmented Lagrangian updates ("<laugon>" - see section 3.3.1 of User's Manual).

Cheers,

Jeff

[smithjh]

Jeff,

To clarify, should the Ogden material with a relatively low bulk modulus work fine for an elastic analysis, a biphasic analysis, or both?

I ask because my student Joel, the original poster, was having trouble getting a biphasic analysis using the Ogden material with a relatively low bulk modulus to converge. Moreover, the failure for the same exact simulation would occur in different elements for different runs.

Thanks,
Josh

[weiss]

Jeff,

To clarify, should the Ogden material with a relatively low bulk modulus work fine for an elastic analysis, a biphasic analysis, or both?

Hi Josh - for an elastic analysis. I will edit my post above to clarify.

I can't think of any reason that a version of the Ogden material model could not be created to use as the solid phase in a biphasic analysis, but a modified version of the current material would have to be implemented.

I believe that uncoupled constitutive models that use the three-field element formulation are not compatible with the biphasic implementation in FEBio. The element formulations (three-field with displacement, pressure and dilation degrees of freedom, and biphasic with displacement and pressure degrees of freedom) are not compatible. Steve or Gerard, could you please confirm? If this is the case, we should add some logic to FEBio to prevent their use. Even better would be to modify the uncoupled materials so that they could be used either with the three-field element formulation or with a standard displacement-based element formulation. If the latter were possible, I think the materials could be used as the solid phase in a biphasic analysis.

Cheers,

Jeff

[ateshian]

Hi Josh,

Thanks for sharing your results with alternative solid matrix constitutive relations. This is also helping me refine my understanding of what might be happening. Let me answer some of your questions:

'I would like to follow-up on your statement that "the solid matrix of biphasic materials is expected to be compressible because the pores of a porous-permeable matrix must be able to change in volume as fluid enters or leaves the matrix." Is this a general rule for biphasic material or just a requirement because of how FEBio is implemented for biphasic materials?'

Basically, fluid flow in a biphasic or poroelastic material is driven by the fluid pressure gradient, according to Darcy's law. In a rigid porous material (the most extreme example of an incompressible material) there can be fluid flow only as a result of an externally applied pressure gradient. In a compressible deformable porous solid, fluid flow can arise either from an externally applied pressure gradient (e.g., a permeation experiment) or as a result of deforming the solid matrix, since the change in pore volume must be accompanied by fluid flow. In an incompressible porous solid, fluid flow can occur only as a result of an externally applied pressure gradient, not as a result of deforming the solid matrix, since the pores cannot change in volume.

Based on these concepts (and your own findings with the Mooney-Rivlin material), I retract my earlier general statement that an incompressible porous solid is incompatible with a biphasic analysis. They should work fine together and they should produce the expected outcome: Under a prescribed pressure gradient, there will be fluid flow in an incompressible porous solid; under an applied deformation (in the absence of an applied pressure gradient), there will not be fluid flow.

I confirmed that FEBio can handle this situation by running an unconfined stress-relaxation analysis of a cylindrical disk: With a compressible solid (neo-Hookean or Holmes-Mow) there is a stress-relaxation response; with in incompressible solid (either Mooney-Rivlin or Ogden), there is no stress-relaxation. It is only in that context that my statement regarding incompatibility of frameworks is meaningful (i.e., we normally expect a biphasic material to exhibit stress-relaxation, but it wouldn't in this case). A user just needs to be aware of these concepts.

'As an aside, I must admit that we are confused by the consistency of your response with a forum posting by Jeff in which he said that "'incompressible' is a mis-label" for the Ogden material and that it "should be fine for compressible values of the bulk modulus."'

What Jeff meant (and of course, he is welcome to comment further on this) is that there is no such thing as a perfectly incompressible material in a finite element implementation. Incompressibility can only be approximated using various implementation methods. Therefore, all implementations of so-called incompressible materials in FEBio (and other finite element codes) are only nearly-incompressible (thus strictly compressible), with the user deciding how much compressibility should be allowed (by adjusting the bulk modulus). So you can use an Ogden material with a zero value for the bulk modulus if you want. However, you should not expect the response to be what Ogden had intended.

'In reviewing the FEBio Theory Manual, I have determined that my misunderstanding came from confusing "solid phase" or "solid constituent" with "solid matrix." I now understand that the solid constituent is incompressible and the solid matrix must be compressible. I further understand why the incompressible Ogden material model is inconsistent with a biphasic formulation.'

The 'skeleton' of a porous solid in a biphasic formulation is indeed intrinsically incompressible. Therefore, applying a hydrostatic fluid pressure to a biphasic material will cause no change in volume, even if the solid is described by a constitutive relation for a compressible material. The compressibility only arises from the change in pore volume. (This is not the most general definition of porous media theories: It is possible to formulate poroelastic materials where the 'skeleton' is intrinsically compressible *and* the pores can change in volume -- however that's not implemented in FEBio and rarely used in biomechanics).

So, having said all of this, I have to assume that the problem you are having with your analysis is not of a fundamental nature. I will take a closer look at your specific analysis and see if I can suggest ways to make it work better. This will take me a little extra time though.

Best,

Gerard

[ateshian]

Hi Josh,

I have attached a file (InflationContact) which conceptually does what I believe you are trying to solve: Two nested hollow biphasic spheres in contact, with the inner sphere subjected to a fluid pressure. The solid matrix is described with an Ogden material where the bulk modulus has a magnitude comparable to the c1 property. The analysis is stable and converges very rapidly.

Here are a few key points that you should know:

When prescribing a fluid pressure on a biphasic surface, you must apply two boundary conditions:

(1) Prescribe the nodal fluid pressure, e.g.,
<prescribe>
<node id="1" bc="p" lc="3">1</node>
...
</prescribe>

(2) Prescribe a surface load of the following type:
<normal_traction type="nonlinear" traction="effective">
<quad4 id="1" lc="2" scale="1"> 1, 59, 303, 55</quad4>
</normal_traction>

This latter boundary condition is documented in the User's Manual but not yet implemented in PreView. In PreView, use Physics->Add Load->Pressure load to prescribe it, then edit the FEBio file manually to replace the <pressure> tag with <normal_traction type="nonlinear" traction="effective"> (and close the tag by replacing </pressure> with </normal_traction>).

Other things to note:
- When performing biphasic analyses, always bias your mesh to create thinner elements near boundaries and contact interfaces. The ability to create biased meshes in hollow spheres is a new feature, you may need to download the latest version of PreView.
- Check out the way I set up the contact interface. I recommend that you turn on autopenalty and laugon (see the attached file).
- It is not necessary to use a multi-step analysis for your problem, you can use must-points instead (see attached file).
- Since I am working with a development version of FEBio, you may have to delete the tag <Module type="poro"/> near the beginning of the file to get your analysis to work.

I hope this example helps. I assume this is the type of analysis you wanted to perform, but please clarify if you need something different.

Best,

Gerard

[smithjh]

Gerard,

I will carefully review your posted simulation. Thank you very much.

Josh

[maas]

Hi all,

I wanted to make a comment on FEBio's implementation of incompressibility. The three-field element formulation that is used in FEBio to enforce incompressibility is optional as of version 1.3.0, although it is always on by default for uncoupled materials. The user can decide not to use the three-field element by setting the following tag in the Control section of the input file.

<use_three_field_hex>0</use_three_field_hex>

However, this only applies to purely elastic simulations. Our biphasic implementation does not use the three-field element formulation, in other words you cannot enforce incompressibility for the solid component of biphasic simulations (unless by increasing the bulk-modulus, but that causes other problems). This has an important consequence: Since the uncoupled materials in FEBio (e.g. Ogden) are designed to be used in a (nearly-)incompressible regime, it may not make much sense to use them in a biphasic analysis. Ideally, we would need to implement alternative compressible formulations of these materials for these types of analyses. Hope this helps.

Cheers,

Steve.

[weiss]

Thanks for clarifying that, Steve!

[smithjh]

Steve,

In conclusion, while FEBio does not prevent the use of incompressible materials in a biphasic analysis (such as what Gerard posted early in this thread), one should not expect physically realistic responses? And, if I am getting good comparisons against analytic solutions for infinitesimal analyses (which is yet to be seen), I would be just lucking out?

What is the likelihood of incorporation of a compressible Ogden material into FEBio that would allow for use in a biphasic analyses?

Thanks,
Josh

[maas]

Hi Josh,

Indeed, the FEBio result may not be physically very realistic. However, that does not imply that good comparisons with analytical solutions cannot be obtained. After all, your analytical result may also be physically unrealistic .

Biphasic analysis is not my specialty and I am not aware of a suitable Ogden formulation. But if you (or anyone else reading this post) know of one, we would be more than happy to add it to FEBio.

Cheers,

Steve.

[smithjh]

Steve,

In Ogden, Non-Linear Elastic Deformations, Dover, 1997, he considers a strain energy function for compressible materials of the form (Eq. 4.4.1 on p. 222)

W = 0.5 mu (lambda_1^2 + lambda_2^2 + lambda_3^2 - 3 - 2 ln J) + 0.5 mu' (J - 1)^2

where mu and mu' are non-zero constants. I would generalize this expression to allow for arbirary exponents, as you do with your incompressible Ogden material.

Thanks,
Josh