" /> FE method for different numbers of elements per mesh.

Figure 9. Deformed shapes for both standard integration (left) and cases, with in-plane shear stress shown. The undeformed mesh is shown in black wireframe.

5.3. Incompressible Mouse Cornea Subjected to Pressure

In the final example, the incompressible behavior of a three-dimensional mouse cornea subjected to pressure loading is studied by means of standard and methods. Biological tissues are prime examples of nearly incompressible media, and the current approach has been applied with good results in [30] [31] . An intraocular pressure of 1.73 kPa is applied to the posterior side of the cornea. The cornea displacement but not rotation is restricted along the center of the cornea edge. While some models of the human cornea use more restrictive boundary conditions on the edges, the mouse cornea is thinner at the edge than at the center, and hence the approximation used here is appropriate. Four different meshes each having 2, 4, 6, and 8 layers of elements through the cornea thickness consisting of respectively 1360, 2720, 4080, and 5440 hexahedral elements are considered to assess convergence.

The mouse cornea is a nearly incompressible soft biological tissue consisting of five layers of which the stroma is the thickest and stiffest. For simplicity, only the stiffness of the stroma is considered in the constitutive model development. The stroma is composed of oriented and dispersed collagen fibrils embedded in nearly incompressible matrix. The material model used is an anisotropic hyperelastic model adapted from Pandolfi and Holzapfel [32] for the human cornea. The model follows the decomposition of the strain energy function into an isochoric component and volumetric components of oriented fibrils distributed in Neo-Hookean matrix:


in which U is the penalty function enforcing incompressibility constraints of the cornea defined as


where is a positive penalty parameter. The isochoric component may be written






for a given unit orientation vector. Here is the shear modulus analogue for the matrix and and are materials parameters that define the stiffening effect of the fibrils. The values of the material parameters used in the simulation are given in Table 5.

The unit vector is the mean orientation of collagen fibrils in the reference configuration. For the human cornea, fibrils are characterized by two mean orientations [32] . However, based on x-ray images [33] , one mean orientation of collagen fibrils is assumed for the mouse cornea and accordingly the quantities involving are removed from the constitutive equations presented in [32] . The mean orientation is assumed to be horizontal (nasal-temporal) near the center, and circumferential near the edge. There is a transition region in between, where the orientation is linearly interpolated. The dispersion parameter describes the ratio of anisotropic fibrils to isotropically distributed ones in the cornea. For the mouse cornea, in the region around the cornea center, we assume that approximately 80% of collagen fibrils are oriented in horizontal direction and 90% have cir-

cumferential orientation around the cornea edge. Again, in the transition region, a linear interpolation of the two directions is assumed.

Spheres with different radii of curvature for outer and inner surfaces of the cornea result in varying thickness throughout. For the mouse cornea model, a thickness of 0.3 mm at the apex and 0.1 mm at the edge was assumed. An initial cornea height of 0.84 mm and the distance from cornea center to the inner and outer edges of respectively 1.33 mm and 1.39 mm were used.

The finite element simulation is performed using the standard 8-node hexahedral elements and the proposed approach, with J averaged in the reference configuration. The convergence is studied by increasing the number of layers of elements through the thickness of the mesh. The quantity of interest is the displacement in the vertical direction at the corneal apex.

The results of vertical displacement at the apex versus the number of element layers through corneal thickness are shown in Figure 10. The results obtained with standard 8-node hexahedral elements and with 8-node elements are compared. Using the method, all the displacement values reach a converged solution with

Figure 10. Vertical Displacement at the cornea apex versus number of layers of elements through the corneal thickness with standard and elements obtained from nonlinear and anisotropic FE simulations of mouse cornea subjected to intraocular pressure.

fewer layers, and the standard elements appear to exhibit some locking. Figure 11 also illustrates the deformed shape of the cornea using the standard and FE simulation. This example demonstrates that the proposed method functions well even with anisotropic materials.

6. Conclusions

In summary, we have developed a method for finite deformation nearly incompressible materials, based on an integral average of the volumetric part of the deformation gradient. There is an advantage in such formulations over reduced-order quadrature formulation in that the deformation can be tracked at fewer integration points, saving memory and perhaps, modestly, computation time. Though not considered here, the advantages are greater in history-dependent materials and multiphysics applications. While integral-averaged methods have been proposed before, there appears to be little in the literature linking a given choice of to the appropriate matrix in these cases, especially for integral formulations. The consistently derived stiffness also ap- pears to be lacking in most of the literature for integral-averaged formulations in standard finite elements, though it has been derived reduced-order quadrature formulations.

We have examined two possible choices for the integral averaging of the Jacobian, over the reference configuration and over the current configuration. While there may be some justifications for the current configuration as more natural, the formulation is more complex. Since the volume change in nearly incompressible materials

Table 5. Material constants assumed for anisotropic and nonlinear FE simulation of mouse cornea.

Figure 11. Vertical displacement mapping of a mouse cornea subjected to pressure obtained from standard and FE simulation (the top and bottom figures respectively) are shown in cross section. The wireframe on the bottom is the undeformed mesh and vertical displacement colored by magnitude is shown on top. The cornea displacement but not rotation is restricted along the center of the cornea edge.

is small, one does not expect there to be a great difference in the results. This observation is confirmed by numerical examples.

Since the formulation is general, it can be applied to other choices of, including reduced integration or logarithmic variations. The formulation is not limited to isotropic materials and can be extended to inelastic materials with some work. Overall, numerical results are similar to existing methods.


The authors gratefully acknowledge the support of the U.S. National Institutes of Health grant 1R21EY020946- 01.

Cite this paper

Craig D.Foster,Talisa MohammadNejad, (2015) Trilinear Hexahedra with Integral-Averaged Volumes for Nearly Incompressible Nonlinear Deformation. Engineering,07,765-788. doi: 10.4236/eng.2015.711067


  1. 1. Hughes, T.J.R. (1987) The Finite Element Method. Prentice-Hall, Upper Saddle River.

  2. 2. Flanagan, D.P. and Belytschko, T. (1981) A Uniform Strain Hexahedron and Quadrilateral with Orthogonal Hourglass Control. International Journal for Numerical Methods in Engineering, 17, 679-706.

  3. 3. Belytschko, T., Ong, J.S.J., Liu, W.K. and Kennedy, J.M. (1984) Hourglass Control in Linear and Nonlinear Problems. Computer Methods in Applied Mechanics and Engineering, 43, 251-276.

  4. 4. Reese, S., Kussner, M. and Reddy, B.D. (1999) A New Stabilization Technique for Finite Elements in Non-Linear Elasticity. International Journal for Numerical Methods in Engineering, 44, 1617-1652.<1617::AID-NME557>3.0.CO;2-X

  5. 5. Reese, S., Wriggers, P. and Reddy, B.D. (2000) A New Locking-Free Brick Element Technique for Large Deformation Problems in Elasticity. Computers & Structures, 75, 291-304.

  6. 6. Sussman, T. and Bathe, K.J. (1987) A Finite Element Formulation for Nonlinear Incompressible Elastic and Inelastic Analysis. Computers & Structures, 26, 357-409.

  7. 7. Chapelle, D. and Bathe, K.J. (1993) The Inf-Sup Test. Computers & Structures, 47, 537-545.

  8. 8. Kasper, E.P. and Taylor, R.L. (2000) A Mixed-Enhanced Strain Method Part I: Geometrically Linear Problems. Computers & Structures, 75, 237-250.

  9. 9. Kasper, E.P. and Taylor, R.L. (2000) A Mixed-Enhanced Strain Method Part II: Geometrically Nonlinear Problems. Computers & Structures, 75, 251-260.

  10. 10. Simo, J.C., Taylor, R.L. and Pister, K.S. (1985) Variational and Projection Methods for the Volume Constraint in Finite Deformation Elasto-Plasticity. Computer Methods in Applied Mechanics and Engineering, 51, 177-208.

  11. 11. Simo, J.C. and Rifai, M.S. (1990) A Class of Mixed Assumed Strain Methods and the Method of Incompatible Modes. International Journal for Numerical Methods in Engineering, 29, 1595-1638.

  12. 12. Masud, A. and Xia, K.M. (2005) A Stabilized Mixed Finite Element Method for Nearly Incompressible Elasticity. Journal of Applied Mechanics-Transactions of the ASME, 72, 711-720.

  13. 13. Masud, A. and Xia, K.M. (2006) A Variational Multiscale Method for Inelasticity: Application to Superelasticity in Shape Memory Alloys. Computer Methods in Applied Mechanics and Engineering, 195, 4512-4531.

  14. 14. Xia, K.M. and Masud, A. (2009) A Stabilized Finite Element Formulation for Finite Deformation Elastoplasticity in Geomechanics. Computers and Geotechnics, 36, 396-405.

  15. 15. Ramesh, B. and Maniatty, A.M. (2005) Stabilized Finite Element Formulation for Elastic-Plastic Finite Deformations. Computer Methods in Applied Mechanics and Engineering, 194, 775-800.

  16. 16. Hughes, T.J.R. (1980) Generalization of Selective Integration Procedures to Anisotropic and Nonlinear Media. International Journal for Numerical Methods in Engineering, 15, 1413-1418.

  17. 17. Brezzi, F. and Fortin, M. (1991) Mixed and Hybrid Finite Element Methods. Springer, Berlin, Heidelerg and New York.

  18. 18. Szabo, B. and Babuska, I. (1991) Finite Element Analysis. Wiley, New York.

  19. 19. Moran, B., Ortiz, M. and Shih, C.F. (1990) Formulation of Implicit Finite-Element Methods for Multiplicative Finite Deformation Plasticity. International Journal for Numerical Methods in Engineering, 29, 483-514.

  20. 20. de Souza Neto, E.A., Peric, D., Dutko, M. and Owen, D.R.J. (1996) Design of Simple Low Order Finite Elements for Large Strain Analysis of Nearly Incompressible Solids. International Journal of Solids and Structures, 33, 3277-3296.

  21. 21. de Souza Neto, E.A., Peric, D. and Owen, D.R.J. (2008) Computational Methods for Plasticity. Wiley, New York.

  22. 22. Nagtegaal, J.C., Parks, D.M. and Rice, J.R. (1974) On Numerically Accurate Finite Element Solutions in the Fully Plastic Range. Computer Methods in Applied Mechanics and Engineering, 4, 153-177.

  23. 23. Mathisen, K.M., Okstad, K.M., Kvamsdal, T. and Raknes, S.B. (2011) Isogeometric Analysis of Finite Deformation Nearly Incompressible Solids. Journal of Structural Mechanics, 44, 260-278.

  24. 24. de Souza Neto, E.A., Pires, F.M.A. and Owen, D.R.J. (2005) F-Bar-Based Linear Triangles and Tetrahedral for Finite Strain Analysis of Nearly Incompressible Solids. Part I: Formulation and Benchmarking. International Journal for Numerical Methods in Engineering, 62, 353-383.

  25. 25. Elguedj, T., Bazilevs, Y., Calo, V.M. and Hughes, T.J.R. (2008) F-Bar Projection Method for Finite Deformation Elasticity and Plasticity Using NURBS Based Isogeometric Analysis. International Journal of Material Forming, 1, 1091- 1094.

  26. 26. Belytschko, T., Tsay, C.S. and Liu, W.K. (1981) A Stabilization Matrix for the Bilinear Mindlin Plate Element. Computer Methods in Applied Mechanics and Engineering, 29, 313-327.

  27. 27. Doll, S., Schweizerhof, K., Hauptmann, R. and Freischläger, C. (2000) On Volumetric Locking of Low-Order Solid and Solid-Shell Elements for Finite Elastoviscoplastic Deformations and Selective Reduced Integration. Engineering Computations, 17, 874-902.

  28. 28. Belytschko, T., Liu, W.K. and Moran, B. (2000) Nonlinear Finite Elements for Continua and Structures. John Wiley & Sons, New York.

  29. 29. Brink, U. and Stein, E. (1996) On Some Mixed Finite Element Methods for Incompressible and Nearly Incompressible Finite Elasticity. Computational Mechanics, 19, 105-119.

  30. 30. Rhee, J., Nejad, T.M., Comets, O., Flannery, S., Gulsoy, E.B., Iannaccone, P. and Foster, C. (2015) Promoting Convergence: The Phi Spiral in Abduction of Mouse Corneal Behaviors. Complexity, 20, 22-38.

  31. 31. Nejad, T.M., Iannaccone, S., Rutherford, W., Iannaccone, P.M. and Foster, C.D. (2015) Mechanics and Spiral Formation in the Rat Cornea. Biomechanics and Modeling in Mechanobiology, 14, 107-122.

  32. 32. Pandolfi, A. and Holzapfel, G.A. (2008) Three-Dimensional Modeling and Computational Analysis of the Human Cornea Considering Distributed Collagen Fibril Orientations. Journal of Biomechanical Engineering-Transactions of the ASME, 130, Article ID: 061006.

  33. 33. Sheppard, J., Hayes, S., Boote, C., Votruba, M. and Meek, K.M. (2010) Changes in Corneal Collagen Architecture during Mouse Postnatal Development. Investigative Ophthalmology & Visual Science, 51, 2936-2942.

Appendix 1: Stiffness Derivation

As mentioned previously, the stiffness matrix may be derived using a directional derivative or a pseudo-time derivative. Here, we take the latter approach. We examine the nodal subvector of the element internal force vector


where, as mentioned previously, and. For preliminaries, note that





. (72)

For convenience, for the remainder of the derivation, we will move the nodal subscript of the and tensors to superscripts, e.g.. Next,



. (75)

The time derivative of the nodal subvector of the element internal force vector, then, can be written

. (76)

Examining the quantity inside the integral

. (77)

We examine the second term first. Recall that is a function of rather than. Hence,


rather than. Continuing





. (83)

Taking the second and third quantities one at a time




. (87)

The third quantity becomes




. (91)

Next we evaluate




. (95)

Adding up all these quantities, we find that


and hence the nodal submatrices of the element stiffness matrix may be written

. (97)

The first term is the common “material stiffness” with the modified matrix, and the second term the standard geometric stiffness. The remaining terms arise from the modifications to the derivatives introduced by modifying the matrix. The final term, of course, depends on the form of. For the Jacobian averaged over the reference configuration






. (103)


. (104)

The process for the Jacobian averaged in the current configuration is the same, though the calculations are more tedious.





. (109)

Hence, the derivative we seek must be

. (110)


*Corresponding author.

Journal Menu >>