# Average of the elastic properties of bio-heterogeneous structures

David Guinovart-Sanjuan, Reinaldo Rodriguez-Ramos, Raul Guinovart-Diaz,

Julian Bravo-Castillero, A. Conci, Frédéric Lebon, Serge Dumont

► To cite this version:

David Guinovart-Sanjuan, Reinaldo Rodriguez-Ramos, Raul Guinovart-Diaz, Julian Bravo-Castillero, A. Conci, et al.. Average of the elastic properties of bio-heterogeneous structures. 23rd ABCM International Congress of Mechanical Engineering - COBEM2015, Dec 2015, Rio de Janeiro, Brazil. pp.CD ROM. hal-01247923

HAL Id: hal-01247923

https://hal.science/hal-01247923

Submitted on 16 Mar 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
AVERAGE OF THE ELASTIC PROPERTIES OF BIO-HETEROGENEOUS
STRUCTURES

David Guinovart Sanjuán

Reinaldo Rodríguez Ramos

Raúl Guinovart Díaz

Julián Bravo Castillero

Facultad de Matemática y Computación, Universidad de La Habana, San Lázaro y L, Vedado, La Habana. CP 10400. Cuba.
davidgs@matcom.uh.cu; reinaldo@matcom.uh.cu; guino@matcom.uh.cu; jbravo@matcom.uh.cu

Andrei Martínez-Finkelshtein

University of Almería and Instituto Carlos I de Física Teórica y Computacional, Granada University, Spain
andrei@ual.es

Aura Conci

Instituto de Computação/UFF, Av. Gal. Milton Tavares de Souza, s/nº, Niterói, CEP: 24210-346 Rio de Janeiro, Brasil.
auraconci2012@gmail.com

Frederic Lebon

Laboratoire de Mécanique et d'Acoustique, Université Aix-Marseille, CNRS, Centrale Marseille, 31 Chemin Joseph-Aiguier, 13402 Marseille Cedex 20, France.
lebon@lma.cnrs-mrs.fr

Serge Dumont

Laboratoire Amiénois de Mathématiques Fondamentale et Appliquées, UPJV, CNRS UMR 7352, 33, rue Saint Leu, 80 039 Amiens Cedex, France.
serge.dumont@u-picardie.fr

**Abstract.** Sometimes the analysis of the stress state, crack initiation and some others phenomenon in heterogeneous media is complicated. In order to facilitate the ulterior study of such problems, we apply the asymptotic homogenization method (AHM) to a general laminated shell composite and an equivalent homogeneous elastic problem with effective properties is obtained. The physical and material properties of the bio-composites are studied using the two scale of AHM. The heterogeneous elastic problem is solved. As examples of the application of the AHM the cornea is considered as a laminated shell bio-composite. The cornea is composed with three cell layers: the outer epithelial cells, stromal cells and endothelial cells. Between these layers are extracellular structures called Bowman's and Descemet's membranes. The cornea has a variable thickness. The polar coordinates system is introduced to study the geometry of this media. A function that describes the variation of the thickness is used to generate the model of the cornea. The effective elastic tensor is calculated by AHM and FEM considering two models, one where the thickness is constants along the cornea and other considering thickness variation. These two models by AHM and FEM are compared and the influence of the variation of the thickness in the effective properties of the cornea is shown. The effective stress is determinate solving the homogenized problem for constant and variable thickness.

**Keywords:** Effective elastic properties, Asymptotic Homogenization Method, Bio-composites, Laminated structures, Numerical analysis, Finite Element Methods

# 1. INTRODUCTION

The cornea possesses a unique combination of mechanical stiffness, strength, and optical transparency that enables it to serve as both a protective covering and the primary refractive component of the eye. These properties depend on the composition and orientation of the cornea, and differ from the central cornea to the limbus, (Nguyen and Boyce, 2011). The stress and the strain are related with the deformation and the geometry of the cornea, hence the importance of the study of the elastic properties of the cornea (Holmes et al., 2001), (Garcia-Porta et al., 2014), (Voorhies, 2003). Otherwise, the glaucoma and myopia are the most common problems that the cornea has on these days (Gatinel et al., 2001), (Broman et al., 2007), (González-Méijope et al., 2008). The transplants and surgery have to be very precise proceeding, then
knowledge of the effective properties become an important tool because in a certain kind of eye surgery, the human eyeball is deformed sustainably by the application of an elastic band. To many author talk about the geometry of the cornea and its influence on the mechanical properties, (Donohue et al., 1995), (Pinsky et al., 2005), (Sigat et al., 2004). Using the two scale Asymptotic Homogenization Method and Finite Elements Method, we determinate the effective properties of the cornea and compare the influence of its geometry on the results.

## 2. LINEAR ELASTIC PROBLEM IN THE CORNEA

In (Skacel and Bursa, 2014), (Genest, 2010), (Li and Tighe, 2006), (Pandolfi and Holzapfel, 2008), (Cardoso and Cowin, 2012), the stress-strain relationship in the cornea has an exponential behavior and the analysis of the elastic properties is made considering a non-linear model. Otherwise, beyond a certain stress, the $\sigma - \epsilon$ relationship becomes almost linear (Elsheikh et al., 2011).

In order to reduce the analysis of the elastic theory problem in the cornea, a linear constitutive equation for the stress and the strain is considered, following the models presented in (Asher et al., 2014), (Cabrera-Fernández et al., 2005), (Donohue et al., 1995), (Fisher, 1971), (Anderson et al., 2004).

The cornea is considered a solid periodic structure $\Omega$ bound by the surfaces $\Sigma_1$ and $\Sigma_2$. Fixed a coordinate system $\mathbf{x} = (x_1, x_2, x_3)$, the geometry of the cornea is described by a function $\rho(\mathbf{x})$. For small deformations, the operator relating the strain $\epsilon$ and the stress $\sigma$ in the cornea is the Hooke's law,

$$ \sigma = C(\mathbf{x}, \mathbf{y}) : \epsilon, \qquad (1) $$

where $C$ is the elastic tensor, regular in $\mathbf{x}$ and periodic in $\mathbf{y}$; here $\mathbf{y} = \rho(\mathbf{x})/\epsilon$ is the fast or local variable and $\epsilon$ is a very small parameter which characterizes the periodicity. Following the Einstein's convention of summing over repeated indices, the elastic problem takes the form

$$ \sigma^{ij} \|_{j} + f_{i} = 0 \text{ on } \Omega, \qquad (2) $$

with boundary conditions

$$ u_j = u_j^0 \text{ on } \Sigma_1, \quad \sigma^{ij} \cdot n_j = S^i \text{ on } \Sigma_2, \qquad (3) $$

where $\{\cdot\}_{\|j}$ denotes the covariant derivative with respect to the variable $x_j$, $f_i$ are the components of the body forces, $u_j$ are the components of the displacement, $n_j$ is the normal vector in $\Sigma_2$ and the functions $u_j^0$ and $S^i$ are the values of the displacement and the stress in $\Sigma_1$ and $\Sigma_2$, respectively.

The general expression of the Cauchy's formula is

$$ \epsilon_{mn} = \frac{1}{2} (u_m \|_n - u_n \|_m). \qquad (4) $$

Taking into account the symmetry of the elastic tensor $C^{ijkl}$ and substituting the equation (4) in equation (1), the contravariant components of the stress tensor are

$$ \sigma^{ij} = C^{ijkl} (u_{k,l} - \Gamma_{kl}^{p} u_{p}), \qquad (5) $$

where $\Gamma_{kl}^p$ denotes the Christoffel's symbols of second kind.

Substituting equation (5) into equation (2), the general expression of the linear elastic problem can be obtained:

$$ \left( \frac{\partial k_{,j}}{\epsilon} C^{ijmn}{}_{|k} + C^{ijmn}{}_{,j} + \Gamma_{jk}^{i} C^{kjmn} + \Gamma_{jk}^{j} C^{ikmn} \right) (u_{m,n} - \Gamma_{mn}^{r} u_{r}) + C^{ijmn} (u_{m,nj} - \Gamma_{mn,j}^{r} u_{r} - \Gamma_{mn}^{r} u_{r,j}) + f^{i} = 0 \quad \text{on } \Omega, \qquad (6) $$

with the boundary conditions

$$ u_i = u_i^0 \text{ on } \Sigma_1 \quad (C^{ijmn} (u_{m,n} - \Gamma_{mn}^r u_r)) n_j = S^i \text{ on } \Sigma_2, \qquad (7) $$

where $\{\cdot\}_{|j} = \frac{\partial}{\partial y_j} \{\cdot\}$ denotes the differentiation with respect to the fast or local curvilinear coordinate.

One of the main difficulties in solving problem (6)-(7) resides in the rapidly oscillating coefficients. In this paper we reduce it to an equivalent problem with homogenized coefficients, and thus obtain the effective properties of the cornea, by using the two-scale Asymptotic Homogenization Method (AHM) described next.
### 3. ASYMPTOTIC HOMOGENIZATION METHOD

The AHM is used to determine the homogenized problem and the effective properties of the solution of (6)-(7). The general expression for the two-scale asymptotic expansion is

$$u_m^{(\epsilon)} = v_m + \epsilon \left[ N_{(0)m}^p v_p + N_{(1)m}^{lk} v_{l,k} \right] + o(\epsilon), \quad (8)$$

where the local functions $N_{(0)m}^{lk} \equiv N_m^{lk}(\mathbf{x}, \mathbf{y})$ and $N_{(0)m}^p \equiv N_{(0)m}^p(\mathbf{x}, \mathbf{y})$ are $\mathbf{y}$-periodic, and $\langle N \rangle = 0$, where $\langle \cdot \rangle$ is the average operator, defined by $\langle \cdot \rangle = \frac{1}{V_Y} \int_Y \langle \cdot \rangle \sqrt{g} d\mathbf{y}$. Here $Y$ is the periodic cell, $V_Y$ is its volume, and $g = \det([g_{ij}])$ where $[g_{ij}]$ is the metric tensor.

Substituting the expansion (8) into the heterogeneous problem (6)-(7) we consider the coefficients for $\epsilon^{-1}$, which must vanish:

$$\left( \varrho_{q,j} C^{ijkl} + \varrho_{p,n} C^{ijmn} N_{(1)m|p}^{lk} \varrho_{q,j} \right) \Big|_{q} = 0, \quad (9)$$

$$\left( -\varrho_{l,j} C^{ijmn} \Gamma_{mn}^{p} + \varrho_{l,j} C^{ijmn} \varrho_{t,n} N_{(0)m|t}^{p} \right) \Big|_{l} = 0. \quad (10)$$

Solving the differential equation (9), the values of the local functions $N_{(1)m|p}^{lk}$ can be determined. Equation (10) allows to calculate $N_{(0)m|t}^p$ as a function of $N_{(1)m|p}^{lk}$,

$$\left( \varrho_{l,j} C^{ijmn} \varrho_{t,n} N_{(0)m|t}^{p} \right) \Big|_{l} = - \Gamma_{mn}^{p} \left( \varrho_{q,s} C^{ijkls} N_{(1)k|q}^{mn} \varrho_{l,j} \right) \Big|_{l}. \quad (11)$$

Applying the average operator to the coefficient of $\epsilon^0$ and grouping the coefficients of $v_{k,jl}$, $v_{k,l}$ and $v_l$, we get

$$C_e^{ijkl} \equiv \hat{h}^{ijkl} = \langle C^{ijkl} + C^{ijmn} \varrho_{p,n} N_{(1)m|p}^{kl} \rangle, \quad (12)$$

as well as expressions for the homogenized coefficients,

$$\hat{h}^{ilk} = \hat{h}_{,j}^{ijkl} - \Gamma_{mn}^{l} \hat{h}^{ikmn} + \Gamma_{jn}^{i} \hat{h}^{njlk} + \Gamma_{jn}^{j} \hat{h}^{inlk}, \quad (13)$$

$$\hat{h}^{il} = - \left( \Gamma_{mn}^{l} \hat{h}^{ijmn} \right)_{,j} - \left( \Gamma_{jk}^{i} \hat{h}^{kjmn} + \Gamma_{jk}^{j} \hat{h}^{ikmn} \right) \Gamma_{mn}^{l}, \quad (14)$$

$$\hat{h}_{*}^{ijl} = -\Gamma_{mn}^{l} \hat{h}^{ijmn}. \quad (15)$$

Using these expressions we write the homogenized problem, consequence of (6)-(7), but now with slowly oscillating coefficients:

$$\hat{h}^{il} v_l + \hat{h}^{ikl} v_{l,k} + \hat{h}^{imln} v_{l,nm} + f^i = 0 \quad \text{on } \Omega, \quad (16)$$

and with boundary conditions

$$v_i = u_i^0 \quad \text{on } \Sigma_1 \left( \hat{h}_{*}^{ijl} v_l + \hat{h}^{ijkl} v_{l,k} \right) n_j = S^i \quad \text{on } \Sigma_2. \quad (17)$$

Now, the effective linear stress has the expression

$$\sigma_e^{ij} = C_e^{ijkl} v_{k,l} + C_e^{ijk} v_k, \quad (18)$$

where $C_e^{ijkl}$ are the effective coefficients (12), $v_k$ is the solution of the homogenized problem (16)-(17) and $C_e^{ijk}$ is

$$C_e^{ijk} = \langle C^{ijmn} N_{(0)m|n}^k - \Gamma_{mn}^k C^{ijmn} \rangle. \quad (19)$$

### 4. COMPOSITION AND GEOMETRY OF THE CORNEA

The cornea is laminate shell composite, with three cell layers: the outer epithelial cells (I), stromal cells (III) and endothelial cells (V). Between these layers are extracellular structures called Bowman's (II) and Descemet's (IV) membranes, (Anderson et al., 2004; Dupps and Wilson, 2006). The layers (I)-(V) of the cornea have central thickness of 14µm, 12µm, 450µm, 10µm and 14µm, respectively. A cross-section of the cornea is shown in Fig. 1.

The spherical coordinate system $(r, \varphi, \theta)$ is used to describe the three-dimensional model of the cornea (Cardoso and Cowin, 2012), shown in Fig. 2a. In order to determine the effective properties of the cornea, we consider it rotationally
<figure><img src="image_2.png" /><figcaption>Figure 1: Cross-section through the cornea</figcaption></figure>

invariant around its symmetry axis. In other words, its elastic properties have no variation in the θ-direction. As a consequence, we can use a two-dimensional model of the cornea, which in turn can be described by means of the polar coordinate system (r, φ).

In our calculations we have used the following values of the cornea parameters: its thickness is approximately $T_1 = 0.5mm$ at its apex and $T_2 = 0.67mm$ at the limbus, (Pandolfi and Holzapfel, 2008; Ruberti et al., 2011); the length of the half cross-section of the cornea is $D_1 = 6mm$, with $D_2 = 2.5mm$ of height. The inner radius in the center of the cornea is $R_1 = 7.56mm$, at the limbus the radius is $R_2 = 7.51mm$ and the angle between the radios is $α = 0.8235$ (see Fig. 2b), (Cabrera-Fernádez et al., 2005).

<figure><img src="image_3.png" /><figcaption>Figure 2: Geometry of the cornea: a) Three dimensional representation of the cornea given by the function (20). b) Geometry of the cross section of the cornea</figcaption></figure>

Taking into account the values of the parameters $R_1, R_2, D_1, D_2, T_1, T_2$ and $α$, we can use the following function that describes the geometry of the cornea:

$$ \rho(r, \varphi) = r + (m \cdot r + n) (\varphi - \pi/2)^2, \qquad (20) $$

with $m = 0.5014, n = -3.8643, r \in [7.56, 8.06]$ and $\varphi \in [\pi/2, \pi/2 + \alpha]$.

<figure><img src="image_4.png" /><figcaption>Figure 3: Comparison of the two-dimensional cornea: thickness varying (continuous line) and constant thickness along the cornea structure (discontinuous line), given by (20) and (21), respectively.</figcaption></figure>

Alternatively, we consider also the constant-thickness model, corresponding to

$$ \rho(r, \varphi) = r. \qquad (21) $$

Using the polar coordinates system both two-dimensional models of the cornea can be compared, see Fig. 3. The positions are given by the equations

$$ x_1 = \rho \cos(\varphi), \quad x_2 = \rho \sin(\varphi), \qquad (22) $$

where $ρ$ is as in (20) and (21).

## 5. EFFECTIVE ELASTIC PROPERTIES OF THE CORNEA

The layers which compose the cornea are isotropic materials; the cornea is oriented un such a way that the fast variation of its properties corresponds to the y-direction, where $y = ρ/ε$ and $ρ$ is given either by (20) or (21). In the
polar coordinate system the metric tensor is $g_{ij} = \text{diag}[1, r^2, 1]$ and the general expression for the isotropic components of $C^{ijkl}$ is

$$C^{ijkl} = \lambda g^{ij} g^{kl} + \mu (g^{ik} g^{jl} + g^{il} g^{jk}), \qquad (23)$$

where $[g^{ij}] = [g_{ij}]^{-1}$.

Considering the polar metric tensor, the nonzero Christoffel's symbols are

$$\Gamma_{22}^{1} = -r \quad \Gamma_{12}^{2} = \Gamma_{21}^{2} = \frac{1}{r}. \qquad (24)$$

The elastic properties for the different layers of the cornea (see Fig. 1) are given in Table 1. These values were taken from (Cabrera-Fernández et al., 2005). The elastic tensor (23) can be determined using the input parameters of Table 1.

<table><caption>Table 1: Mechanical properties of the cornea.</caption><thead><tr><th>Layers</th><th>Young's Modulus</th><th>Poisson's ratio</th></tr></thead><tbody><tr><td>Epithelial cells (I)</td><td>0.622</td><td>0.49</td></tr><tr><td>Bowman's membrane (II)</td><td>0.275</td><td>0.49</td></tr><tr><td>Stromal cells (III)</td><td>0.279</td><td>0.49</td></tr><tr><td>Descemet's membrane (IV)</td><td>0.304</td><td>0.49</td></tr><tr><td>Endothelial cells (V)</td><td>0.331</td><td>0.49</td></tr></tbody></table>

Considering $y = \varrho/\epsilon$ one-dimensional, the local problems (9) and (10) take the following expressions, respectively:

$$\frac{\partial}{\partial y} \left( \varrho_{,j} C^{ijkl} + \varrho_{,n} C^{ijmn} \frac{\partial N^{lk}_{(1)m}}{\partial y} \varrho_{,j} \right) = 0, \qquad (25)$$

$$\frac{\partial}{\partial y} \left( -\varrho_{,j} C^{ijmn} \Gamma_{mn}^{p} + \varrho_{,j} C^{ijmn} \varrho_{,n} \frac{\partial N^{p}_{(0)m}}{\partial y} \right) = 0. \qquad (26)$$

In order to obtain the effective coefficients, the value of the local function $\partial N^{lk}_{(1)m}/\partial y$ is determined by integrating the local problem (25),

$$\frac{\partial N^{lk}_{(1)m}}{\partial y} = (A^{ilk} - \varrho_{,p} C^{iplk}) (\varrho_{,n} C^{ijmn} \varrho_{,j})^{-1}, \qquad (27)$$

where the tensor $A^{ilk} \equiv A^{ilk}(\mathbf{x})$. With account of the periodic condition of the local function $N^{lk}_m$, the average $\langle \partial N^{lk}_{(1)m}/\partial y \rangle = 0$, and the function $A^{ilk}$ yields

$$A^{ilk} = \langle \varrho_{,p} C^{iplk} (\varrho_{,n} C^{ijmn} \varrho_{,j})^{-1} \rangle \langle (\varrho_{,n} C^{ijmn} \varrho_{,j})^{-1} \rangle^{-1}. \qquad (28)$$

Finally, functions $\partial N^{lk}_{(1)m}/\partial y$ are obtained by using the values of $\varrho$ and $C^{ijkl}$, given in (20) and (23), respectively, in the equations (27) and (28).

In this paragraph, the numerical method proposed to solve problems (25), (26) is described. Since the problem is composed of two coupled PDE's with singular (discontinuous) rigidity tensor, a $P_1$ (piecewise linear) finite element method (FEM) is used on a weak formulation of (25), (26).

Due to the properties of the finite element discretization, after the computation of an approximation of the solution $N^{p}_{(0)m}$ and $N^{lk}_{(1)m}$ using the FEM method, the integrals (see formula (12) for example) for the calculus of the effective coefficients are replaced by sums over each element contribution.

In what follows we use the super-indices $t_v$ and $t_0$ in order to denote the computed coefficients corresponding to the thickness functions (20) and (21), respectively.

In (Anderson et al., 2004), (Asher et al., 2014), (Pandolfi and Holzapfel, 2008), (Lanchares et al., 2008), (Li and Tighe, 2006) and (Ruberti et al., 2011) the variable thickness along animal or human lens is considered. Otherwise, (Broman et al., 2007), (Last et al., 2012), (Donohue et al., 1995) and (Huang et al., 2002) a constant average thickness in the cornea's structure is analyzed. In this work, we address and compare both cases in order to analyze the influence of the thickness variation on the effective properties.

We also compare the values of the effective coefficients $C^{ijkl}$ obtained by AHM and the finite element method (FEM) in the both models. The effective coefficients for $t_0$, found by AHM and FEM, are shown in Table 2 for the values $r = 7.56$ and $\varphi = \pi/2$.
<table><caption>Table 2: Effective coefficients $C_e^{t_0}$ calculated by AHM and FEM. Difference between $C_e^{t_0}$ and $C_e^{t_v}$.</caption><thead><tr><th>Coefficients</th><th>C<sup>t_0</sup><sub>e</sub> (AHM)</th><th>C<sup>t_0</sup><sub>e</sub> (FEM)</th><th>||C<sup>t_0</sup><sub>e</sub> - C<sup>t_v</sup><sub>e</sub>|| (AHM)</th></tr></thead><tbody><tr><td>C<sup>1111</sup><sub>e</sub></td><td>4.8812</td><td>4.8814</td><td>8 × 10<sup>-4</sup></td></tr><tr><td>C<sup>1122</sup><sub>e</sub></td><td>0.0821</td><td>0.0836</td><td>8 × 10<sup>-8</sup></td></tr><tr><td>C<sup>1133</sup><sub>e</sub></td><td>4.6898</td><td>4.6936</td><td>2 × 10<sup>-6</sup></td></tr><tr><td>C<sup>2222</sup><sub>e</sub></td><td>0.0015</td><td>0.0015</td><td>2 × 10<sup>-9</sup></td></tr><tr><td>C<sup>2233</sup><sub>e</sub></td><td>0.0821</td><td>0.0821</td><td>4 × 10<sup>-8</sup></td></tr><tr><td>C<sup>3333</sup><sub>e</sub></td><td>4.8886</td><td>4.8888</td><td>3 × 10<sup>-15</sup></td></tr><tr><td>C<sup>1112</sup><sub>e</sub></td><td>0</td><td>0</td><td>1 × 10<sup>-8</sup></td></tr><tr><td>C<sup>2212</sup><sub>e</sub></td><td>0</td><td>0</td><td>3 × 10<sup>-7</sup></td></tr><tr><td>C<sup>3312</sup><sub>e</sub></td><td>0</td><td>0</td><td>1 × 10<sup>-5</sup></td></tr><tr><td>C<sup>2323</sup><sub>e</sub></td><td>0.0017</td><td>0.0017</td><td>2 × 10<sup>-8</sup></td></tr><tr><td>C<sup>1323</sup><sub>e</sub></td><td>0</td><td>0</td><td>6 × 10<sup>-6</sup></td></tr><tr><td>C<sup>1313</sup><sub>e</sub></td><td>0.0957</td><td>0.0957</td><td>1 × 10<sup>-6</sup></td></tr><tr><td>C<sup>1212</sup><sub>e</sub></td><td>0.0017</td><td>0.0017</td><td>8 × 10<sup>-8</sup></td></tr></tbody></table>

In the case of the varying thickness (denoted by the superscript $t_v$), the coefficients do depend on $r$ and $\varphi$ (see Tables 3 and 4). Thus, in Table 2 we present also the distance between $C_e^{t_0}$ and $C_e^{t_v}$.

We present also the values of the components $C_e^{1111}$ (Table 3) and $C_e^{1112}$ (Table 4) of the elastic tensor $C_e^{t_v}$ for the case of the varying thickness.

<table><caption>Table 3: Components $C_e^{1111}$ of the elastic tensor $C_e^{t_v}$ for different values of $r$ and $\varphi$.</caption><thead><tr><th rowspan="2">r ↓</th><th colspan="2">1.5708</th><th colspan="2">1.9002</th><th colspan="2">2.3943</th></tr><tr><th>AHM</th><th>FEM</th><th>AHM</th><th>FEM</th><th>AHM</th><th>FEM</th></tr></thead><tbody><tr><td>7.56</td><td>4.88123</td><td>4.88133</td><td>4.88098</td><td>4.88132</td><td>4.88083</td><td>4.88129</td></tr><tr><td>7.76</td><td>4.88123</td><td>4.88133</td><td>4.88132</td><td>4.88134</td><td>4.88137</td><td>4.88133</td></tr><tr><td>8.06</td><td>4.88123</td><td>4.88133</td><td>4.88182</td><td>4.88127</td><td>4.88211</td><td>4.88221</td></tr></tbody></table>

<table><caption>Table 4: Components $C_e^{1112}$ of the elastic tensor $C_e^{t_v}$ for different values of $r$ and $\varphi$.</caption><thead><tr><th rowspan="2">r ↓</th><th colspan="2">1.5708</th><th colspan="2">1.9002</th><th colspan="2">2.3943</th></tr><tr><th>AHM</th><th>FEM</th><th>AHM</th><th>FEM</th><th>AHM</th><th>FEM</th></tr></thead><tbody><tr><td>7.56</td><td>0</td><td>0</td><td>0.0243 × 10<sup>-4</sup></td><td>0.0272 × 10<sup>-4</sup></td><td>0.0479 × 10<sup>-4</sup></td><td>0.0535 × 10<sup>-4</sup></td></tr><tr><td>7.76</td><td>0</td><td>0</td><td>-0.0088 × 10<sup>-4</sup></td><td>-0.0098 × 10<sup>-4</sup></td><td>-0.0172 × 10<sup>-4</sup></td><td>-0.0193 × 10<sup>-4</sup></td></tr><tr><td>8.06</td><td>0</td><td>0</td><td>-0.0584 × 10<sup>-4</sup></td><td>-0.0653 × 10<sup>-4</sup></td><td>-0.1148 × 10<sup>-4</sup></td><td>-0.1285 × 10<sup>-4</sup></td></tr></tbody></table>

We can appreciate a very good coincidence between the effective properties calculated by AHM and FEM. It should be taken into account that in the present situation even a small difference in the effective coefficients (for varying and constant thicknesses), as in the third column of Table 2, renders important difference in the anisotropy of the material. Namely, for the constant thickness with isotropic layers the homogenized media exhibits transversely isotropic behavior, but in the case with varying thickness the material has a monoclinic symmetry.

A important comparison between the calculation of the effective coefficients, obtained by AHM and FEM, is the time of the computational operations. The values of Table 2 have been reproduced using AHM in an average time of 0.051 seconds, however, the average time reported by FEM was 1.6 seconds.

## 6. SOLVING THE HOMOGENIZED PROBLEM OF THE CORNEA

An external uniform force is applied on a two-dimensional cornea, simulating the intra-ocular pressure (IOP). This force $f$ occurs in the $r$-direction and has contravariant components $(\cos \varphi, 0, 0)$ with respect to the basis of the tangent vectors $g_i$. This force creates an internal displacement $u|_{\Sigma_1} = (0, 0, 0)$ and the traction $S|_{\Sigma_2}$ that has the same contravariants components of the force $f$.

The influence of the thickness variation in the stress tensor is analyzed. We consider a particular case of the elastic problem. The variation of the properties is in the $x_1$ direction and the only non zero displacement is $u_1$. Under these assumptions, the homogenized problem (16)-(17) is

$$ \hat{h}^{11}v_1 + \hat{h}^{111}v_{1,1} + \hat{h}^{1111}v_{1,11} + \cos \varphi = 0 \quad \text{on} \ r \in (7.56, 8.06), \qquad (29) $$
with the boundary conditions

$$v_1 = 0 \text{ on } r = 7.56 \left( \hat{h}_*^{111} v_1 + \hat{h}^{1111} v_{1,1} \right) n_1 = \cos \varphi \text{ on } r = 8.06. \quad (30)$$

for $\varphi \in [\pi/2, \pi/2 + \alpha]$.

The resulting coefficients $\hat{h}^{1111}$ appear in Table 2; the corresponding coefficients $\hat{h}^{11}$, $\hat{h}^{111}$ and $\hat{h}_*^{111}$ are determined using equations (13)-(15).

The problem (29)-(30) is solved at different points of the cornea. The effective stress is calculated using the equations (18) and (19) for the angles $\varphi = \pi/2$ (corneal apex) $\varphi = 1.9002$ (a point between the apex and the limbus) and $\varphi = \pi/2 + \alpha$ (limbus of the cornea).

<figure><img src="image_5.png" /><figcaption>Figure 4: Comparison of the behavior of the effective stress $\sigma_e^{11}$ at different points of the cornea with constant and variable thickness along the structure.</figcaption></figure>

The component $\sigma_e^{11}$ of the effective stress is varying along the cornea and has a non-linear behavior. At the center of the cornea ($\varphi = 1.5708$) the component $\sigma_e^{11}$ for variable ($\diamond$) and constant thickness($-$) coincide. However, the effective stress has different behavior for other values of $\varphi$, effect being more apparent at the limbus, ($\varphi = 2.2778$). Notice that the highest stress is located at the center of the cornea and the lowest one corresponds to the limbus.

## 7. CONCLUSIONS

The AHM is a powerful tool to determine the effective properties of the cornea quickly and with high precision. The methodology bring the possibility to aboard problems like the analysis of the bio-mechanics properties of the cornea after an intra-stromal ring segments implantation (Martínez-Finkelshtein et al., 2009), (Martínez-Finkelshtein et al., 2011), or the study of the variation of the corneal thickness in a healthy eye as function of the anterior face of the cornea and inner ocular pressure (IOP).

The difference of the length along the cornea is studied and its influence in the mechanical behavior of the cornea as a laminated shell bio-composite. The analytic expressions of the effective coefficients are very important result of this work. The numerical values are compared with the effective coefficients obtain by FEM, very good coincidence is reported and significant difference between the time of the computational operations.

## ACKNOWLEDGEMENTS

The authors gratefully acknowledge to the project SHICHAN, supported by FSP (Cooperation Scientifique Franco-Cubaine) Cuba 2011-26, PROJET N° 29935XH and to the project Composite Materials from University of Havana. RRR thanks to Department of International Relations - DRI - Universidade Federal Fluminense - UFF, Edital Programa de Apoio a Projetos de Internacionalização da UFF - PIUFF - Edição 2014 for the financial support during his stay at UFF as a visiting professor in 2015. AC gratefully thank the institutional support provided by INCT-MACC (National Institute of Science and Technology in Computer Science Medical Assistance and the Brazilian agency CNPq Project number 302298/2012-6. The fifth author (AMF) was partially supported by MICINN of Spain and by the European Regional Development Fund (ERDF) under grant MTM2011-28952-C02-01, by Junta de Andalucía (Excellence Grant P11-FQM-7276 and the research group FQM-229), and by Campus de Excelencia Internacional del Mar (CEIMAR) of the University of Almería. This work was completed during a visit of AMF to the Department of Mathematics of the Havana University, whose hospitality he wishes to acknowledge.
REFERENCES

Anderson, K., El-Sheikh, A. and Newson, T., 2004. "Application of structural analysis to the mechanical behaviour of the cornea". *The Royal Society*, Vol. 1, pp. 3–15.

Asher, R., Gefen, A., Moisseiev, E. and Varssano, D., 2014. "Etiology of keratoconus: proposed biomechanical pathogenesis". *In Silico Cell and Tissue Science*, Vol. 1:3, pp. 1–9.

Broman, A.T., Congdon, N.G., Bardeen-Roche, K. and Quigley, H.A., 2007. "Influence of corneal structure, corneal responsiveness, and other ocular parameters on tonometric measurement of intraocular pressure". *J Glaucoma*, Vol. 16, pp. 581–588.

Cabrera-Fernádez, D., Naizy, A.M., Kurtz, R., Djotyan, G.P. and Juhasz, T., 2005. "Finite element analysis applied to cornea reshaping". *Journal of Biomedical Optics*, Vol. 10(6), pp. 064018–064018–11.

Cardoso, L. and Cowin, S.C., 2012. "Role of structural anisotropy of biological tissues in poroelastic wave propagation". *Mechanics of Materials*, Vol. 44, pp. 174–188.

Donohue, D.J., Stoyanov, B.J., McCally, R.L. and Farrell, R.A., 1995. "Numerical modeling of the cornea's lamellar structure and birefringence properties". *Optical Society of America*, Vol. 12, pp. 1425–1438.

Dupps, W.J. and Wilson, S.E., 2006. "Biomechanics and wound healing in the cornea". *Experimental Eye Research*, Vol. 83, pp. 709–720.

Elsheikh, A., Kassem, W. and Jones, S.W., 2011. "Strain-rate sensitivity of porcine and ovine corneas". *Acta of Bioengineering and Biomechanics*, Vol. 13, pp. 25–36.

Fisher, R.F., 1971. "The elastic constants of the human lens". *J. Physiol*, Vol. 212, pp. 147–180.

Garcia-Porta, N., Fernandes, P., Queiros, A., Salgado-Borges, J., Parafita-Mato, M. and González-Méijome, J.M., 2014. "Corneal biomechanical properties in different ocular conditions and new measurement techniques". *ISRN Ophthalmology*, Vol. 2014, pp. 1–19.

Gatinel, D., Hoang-Xuan, T. and Azar, D.T., 2001. "Determination of corneal asphericity after myopia surgery with the excimer laser: A mathematical model". *Investigative Ophthalmology & Visual Science*, Vol. 42, pp. 1736–1742.

Genest, R., 2010. "Effect of intraocular pressure on chick eye geometry, finite element modeling, and myopia".

González-Méijope, J.M., Villa-Collar, C., Queirós, A., Jorge, J. and Parafita, M.A., 2008. "Study on the influence of corneal biomechanical properties over the short term in response to corneal refractive therapy for myopia". *Cornea*, Vol. 27, pp. 421–426.

Holmes, D.F., Gilpin, C.J., Baldock, C., Ziese, U., Koster, A.J. and Kadler, K.E., 2001. "Corneal collagen fibril structure in three dimensions: Structural insights into fibril assembly, mechanical properties, and tissue organization". *PNAS*, Vol. 98, pp. 7307–7312.

Huang, D., Tang, M. and Shekhar, R., 2002. "Mathematical model of corneal surface smoothing after laser refractive surgery". *Am J Ophthalmol*, Vol. 135, pp. 267–278.

Lanchores, A., Calvo, B., Cristóbal, J.A. and Doblaré, M., 2008. "Finite element simulation of arcuates for astigmatism correction". *Journal of Biomechanics*, Vol. 41, pp. 797–805.

Last, J.A., Thomasy, S.M., Croasdale, C.R., Russell, P. and Murphy, C.J., 2012. "Compliance profile of the human cornea as measured by atomic force microscopy". *Micron*, Vol. 43, pp. 1293–1298.

Li, L. and Tighe, B., 2006. "Nonlinear analysis of static axisymmetric deformation of the human cornea". *Computational Materials Science*, Vol. 38, pp. 618–624.

Martínez-Finkelshtein, A., Delgado, A.M., Castro, G.M., Zarzo, A. and Alió, J.L., 2009. "Comparative analysis of some modal reconstruction methods of the shape of the cornea from the corneal elevation data". *Investigative Ophthalmology and Visual Science*, Vol. 50, No. 12, pp. 5639–5645.

Martínez-Finkelshtein, A., Ramos-López, D., Castro-Luna, G.M. and Alió, J.L., 2011. "Adaptive cornea modeling from keratometric data". *Investigative Ophthalmology and Visual Science*, Vol. 52, No. 8, pp. 4963–4970.

Nguyen, T.D. and Boyce, B.L., 2011. "An inverse finite element method for determining the anisotropic properties of the cornea". *Biomech Model Mechanobiol*, Vol. 10, pp. 323–337.

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*, Vol. 130, pp. 061006–1–061006–12.

Pinsky, P.M., van der Heide, D. and Chernyak, D., 2005. "Computational modeling of mechanical anisotropy in the cornea and sclera". *J Cataract Refract Surg*, Vol. 31, pp. 136–145.

Ruberti, J.W., Sinha-Roy, A. and Roberts, C.J., 2011. "Corneal biomechanics and biomaterials". *Annual Review of Biomedical Engineering*, Vol. 13, pp. 269–295.

Sigat, I.A., Flanagan, J.G., Tertinegg, I. and Ethier, C.R., 2004. "Finite element modeling of optic nerve head biomechanics". *Investigative Ophthalmology & Visual Science*, Vol. 45, pp. 4378–4387.

Skacel, P. and Bursa, J., 2014. "Comparison of constitutive models of arterial layers with distributed collagen fibre orientations". *Acta of Bioengineering and Biomechanics*, Vol. 16, pp. 47–58.

Voorhies, K.D., 2003. "Static and dynamic stress/strain properties for human and porcine eyes".