La formulación de hiperelasticidad incompresible en la matriz de rigidez general del método de los elementos finitos

Contexto

Los pasados cinco años he estado bien involucrado en el estudio del modelado constitutivo y numérico de materiales elásticos no-lineales, entre los que se pueden encontrar materiales tipo-goma, materia biológica como tejidos y polímeros naturales, así como una gran gama de polímeros sintéticos que hallan comportamientos elásticos a grandes deformaciones.

Figura. Simulación numérica de la compresión de una pelota de material elástico no-lineal.

La cuestión es que, cuando comencé el estudio de esto, confronté con que la información está dispersa en un sinfín de fuentes y referencias. Por lo que, al menos en mi experiencia y al carecer de un guía experto en la materia, tuve a bien el tener que hacer una revisión pieza por pieza de los elementos que conforman ese gran rompecabezas del modelado constitutivo de los materiales que se mencionan arriba.

Esta entrada será devota a hacer mención, a un alto nivel, sobre cómo es que se formula a través del Elemento Finito uno de los herramentales más socorridos dentro del modelado de materiales elásticos no-lineales: la hiperelásticidad. Con todo, espero que sea útil sólo como revisión al concepto y el lenguaje utilizado en la implementación de este. No pretendo desarrollar código, dado que en la actualidad existe un sinfín de aplicaciones y software, tanto libre y comercial, que tiene ya implementado esto (ver ABAQUS, FEBio, ANSYS).

En realidad, la pretensión detrás de esta entrada es más bien filosófica, en términos de exponer la coherencia lingüística del modelado de elementos finitos hiperelásticos. Se sugiere, como siempre, hacer una revisión a las fuentes bibliográficas apuntadas al final de la entrada.

Hiperelásticidad

Un material en el cual su respuesta a la deformación, i.e., su tensor de esfuerzos, se deriva de una función de energía de deformación y, además, su deformación ocurre sin ningún cambio apreciable en el volumen al mismo tiempo que posee la capacidad de ser reversible termodinámicamente, se le conoce como material hiperelástico incompresible. Esto, en un sentido más general implica que dicha respuesta sigue siendo elástica, incluso a grandes deformaciones, por lo que se presenta la consecuente dificultad cuando se modela en variantes de deformación plana y axisimetría, pues el tensor de esfuerzos no se determina sólo a través de la deformación. Por lo anterior, se debe introducir una presión hidrostática para evaluar dicho tensor de esfuerzos, sin afectar, al mismo tiempo, al tensor de deformación. A tal presión se le refiere como multiplicador de Lagrange.

La función de energía de deformación para materiales hiperelásticos propuesta por Mooney y Rivlin, $U$, posee la siguiente forma:

\begin{equation}
\bar{U}=U(I_{1}, I_{2})+p(I_{3}-1),
\end{equation}

(1)

donde $p$ es el multiplicador de Lagrange, $I_{1}$, $I_{2}$ y $I_{3}$ son los invariantes de deformación del tensor de deformación de Green. Aunque existen diversas formas de $\bar{U}$, de complejidad y sofisticado, en este ensayo se menciona al modelo de Mooney-Rivlin sólo para efectos ilustrativos y se recomienda consultar el libro de Bower (2009) para considerar un mayor alcance y conocimiento de las mismas.

$\bar{U}$ se usa para el desarrollo de la matriz de rigidez de un elemento y armar una matriz de rigidez global de un sistema tomando en cuenta el comportamiento de un material hiperelástico. Hasta aquí, de manera muy resumida, lo que es un material hiperelástico y

Formulación de hiperelásticidad en la Matriz de Rigidez Global de los Elementos Finitos

Por definición, para un material hiperelástico, el segundo tensor de esfuerzo de Piola-Kirchhoff, $S_{ij}$, se obtieve al derivar la Ec. (1) con respecto a componentes del tensor de Green-Lagrange, $E_{ij}$, como sigue:

\begin{equation}
S_{ij}=\bigg[\frac{\partial U(E_{kl})}{\partial E_{ij}}\bigg]+\bigg[\frac{\partial I_3(E_{kl})}{\partial E_{ij}}\bigg] .
\end{equation}

(2)

Se utilizan las Ecs. (1) y (2) para la formulación de un elemento finito, arbitrario, de cuerpo continuo y perteneciente a una configuración requerida, de acuerdo con el Principio de los Desplazamientos Virtuales [PDV]:

\begin{equation}
\int_{V_0}^{}S_{ij}\delta E_{ij} dV= \delta U^{e},
\end{equation}

(3)

donde $V_0$ es el volumen del elemento en su configuración inicial, y $\delta U^{e}$ es un trabajo virtual externo.

Sustituyendo (2) en (3), aplicando la expansión de Taylor para una configuración previamente calculada (i.e, $\bar{E_{ij}}$, $\bar{S_{ij}}$, $\bar{p}$), y negando términos no-lineales de las incógnitas, se desarrollan las ecuaciones lineales para desplazamientos incrementales y el multiplicador $p$ como sigue:

\begin{eqnarray}
\int_{V_0}^{}\Delta e_{rs} C_{ijrs}\delta \Delta e_{ij} dV +\Delta p\int_{V_0}^{}\frac{{\partial I_3 (E_{kl})}}{{\partial E_{ij}}}\delta \Delta e_{ij} dV+\int_{V_0}^{} \bar{S_{ij}}\delta \Delta \eta_{ij} dV \nonumber \\
= \delta U^{e} - \int_{V_0}^{} \bar{S_{ij}}\delta \Delta e_{ij} dV ,
\end{eqnarray}

(4)

donde $\Delta e_{ij}$ y $\Delta \eta_{ij}$ son componentes lineales y no-lineales incrementales del tensor de deformación, $\Delta p$ es un incremental del multiplicador de Lagrange, y:

\begin{equation} \label{C}
C_{ijrs} =\bigg[\frac{\partial{U(\bar{E}_{kl})}}{\partial {E_{ij}} \partial{E_{rs}}}\bigg] .
\end{equation}

(5)

Ya que la Ec. (4) incluye una variable extra, $\Delta p$, la condición de incompresibilidad se satisface sobre todo el elemento de manera lineal de acuerdo con lo siguiente:

\begin{equation}
\int_{V_0}^{} \frac{\partial{I_3(\bar{E}_{kl})}}{\partial{E_{ij}}} \Delta e_{ij} dV= -\int_{V_0}^{}\bigg[I_3(\bar{E}_{kl}-1)\bigg]dV .
\end{equation}

(6)

Es importante enfatizar que en las anteriores ecuaciones, la presión hidrostática se ha asumido como constante sobre todo el elemento, por que se espera que una deformación constante sea acorde con el problema de estudio de un material hiperelástico.

Al rescribir las Ecuaciones. (4) y (6) en forma matricial simple, fácilmente se puede visualizar la matriz de rigidez característica del elemento en su forma simétrica:

\begin{equation} \label{stiffness}
\quad
\begin{bmatrix}
k_{L}+k_{NL} & k_{I} \\
{k^{T}_{I}} & 0
\end{bmatrix}
\begin{Bmatrix}
\Delta u \\
\Delta p
\end{Bmatrix}
=\begin{Bmatrix}
r-f \\
i
\end{Bmatrix},
\end{equation}

(7)

donde:

•  $[k_L] = \int_{V_0}^{} [B_{L}]^T [C][B_{L}]dV$,
•  $[k_{NL}] = \int_{V_0}^{} [B_{NL}]^T [\bar{S}][B_{NL}]dV$,
•  $\{k_I\} = \int_{V_0}^{} [B_{L}]^T\big\{{\frac{{I_{3}(\bar{E}_{kl})}}{\partial{E_{ij}}}}\big\}dV$,
•  $\{f\} = \int_{V_0}^{} [B_{L}]^T \{\ \bar{S}\}dV$,
•  $i = -\int_{V_0}^{} [I_{3}(E_{kl})-1]dV$.

$\{ \Delta u\}$ y $\{ r\}$ son las columnas de los desplazamientos nodales y las fuerzas, $\{ \Delta e\} = [B_{L}]\{ \Delta u\}$ y $[B_{NL}]$ y $[\bar{S}]$ se obtienen del Método de la Matriz de Rígidez (MMR) (Bathe, 2006).

Al ensamblar las matrices de una estructura, siguiendo el procedimiento estándar del Método de los Elementos Finitos, se obtendría:

\begin{equation}
\quad
\begin{bmatrix}
K_{L}+K_{NL} & K_{I} \\
{K^{T}_{I}} & 0
\end{bmatrix}
\begin{Bmatrix}
\Delta U \\
\Delta P
\end{Bmatrix}
=\begin{Bmatrix}
R-F \\
I
\end{Bmatrix} ,
\end{equation}

(8)

que corresponde con la matriz de rigidez global de un sistema dado.

Comentarios finales

Un material hiperelástico incompresible debe cumplir dos condiciones: (i) su respuesta mecánica a la deformación puede derivarse de una función de energía de deformación, y (ii) al verse deformado, su volumen se conserva. A partir de ahí, se tiene que el segundo tensor de Piola-Kirchhoff se obtiene derivando una función energética respecto al tensor de deformación de Green y al introducir un multiplicador de Lagrange, con lo que subsecuentemente se puede desarrollar el PDV para concluir con la formulación de la matriz de rigidez de un elemento finito dado. Finalmente, aplicando el procedimiento del Método de los Elemento Finitos es posible armar una matriz de rigidez global de un sistema.

Lo anterior permite concluir lo siguiente: la integración de la teoría de hiperelásticidad para materiales no compresibles y el método de los elementos finitos se puede entender, a un alto nivel, sí se precisa como el involucramiento de un potencial energético de deformación para el desarrollo de las matrices de rigidez individuales de un elemento tal como se muestra en la Ecuación (4).

Muchos polímeros, especialmente los elastómeros, pueden ser modelados a través de sistemas como el que se presenta en la Ecuación (8).

Gracias al avance y desarrollo de teorías más sofisticadas de hiperelásticidad y con el avance tecnológico, existen muchos códigos que permiten al usuario implementar rutinas para el modelo de la respuesta de materiales no-lineales como los elastómeros o materiales biológicos, sin embargo, es importante precisar el conocimiento general que se oculta en dichos códigos. Primariamente a un alto nivel y, en la medida de lo posible, de manera detallada.

Bibliografía sugerida

• Bathe, K. (2006). Finite Element Procedures. Prentice Hall.
• Bower, A. (2009). Applied Mechanics of Solids. Taylor & Francis Ed.