scirp.org/file/4-8102473x222.png" /> 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   . 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  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  . However, based on x-ray images  , 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  . 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.
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
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
For convenience, for the remainder of the derivation, we will move the nodal subscript of the and tensors to superscripts, e.g.. Next,
The time derivative of the nodal subvector of the element internal force vector, then, can be written
Examining the quantity inside the integral
We examine the second term first. Recall that is a function of rather than. Hence,
rather than. Continuing
Taking the second and third quantities one at a time
The third quantity becomes
Next we evaluate
Adding up all these quantities, we find that
and hence the nodal submatrices of the element stiffness matrix may be written
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
The process for the Jacobian averaged in the current configuration is the same, though the calculations are more tedious.
Hence, the derivative we seek must be