Caracterización mecánica y modelización hiperelástica de hidrogeles de colágeno en aplicaciones de mecanobiología celular

Proyecto Fin de Carrera Ingeniería de Telecomunicación Formato de Publicación de la Escuela Técnica Superior de Ingeniería Autor: F. Javier Payán Somet Tutor: Juan José Murillo Fuentes Dep. Teoría de la Señal y Comunicaciones Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2013 Trabajo Fin de Máster Máster Universitario en Ingeniería Industrial Caracterización mecánica y modelización hiperelástica de hidrogeles de colágeno en aplicaciones de mecanobiología celular Autor: Alejandro Apolinar Fernández Tutor: José Antonio Sanz Herrera Departamento de Mecánica de Medios Continuos y Teoría de Estructuras Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2020 Trabajo Fin de Máster Máster Universitario en Ingeniería Industrial Caracterización mecánica y modelización hiperelástica de hidrogeles de colágeno en aplicaciones de mecanobiología celular Autor: Alejandro Apolinar Fernández Tutor: José Antonio Sanz Herrera Profesor Titular Departamento de Mecánica de Medios Continuos y Teoría de Estructuras Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2020 A mis padres, Gustavo y Verónica, por su eterno apoyo, sin el cual este trabajo no habría sido posible. Resumen En la consecución de este trabajo se ha estudiado la viabilidad del empleo de tres modelos constitutivos hiperelásticos en la caracterización de hidrogeles de colágeno: modelo Neo-Hookeano, modelo de Holzapfel y modelo de Steinwachs. Estos geles son utilizados en mecanobiología para simular las características de la matriz extracelular, siendo de gran importancia en el marco de la Microscopía de Fuerza de Tracción (TFM, por sus siglas en inglés). En dicho estudio, se ha efectuado la implementación de los modelos en MATLAB, y sirviéndose de las herramientas que el programa ofrece, se han realizado ajustes por mínimos cuadrados de los parámetros gobernantes de dichos modelos a curvas experimentales enviadas por colaboradores de la Universidad KU Leuven. El modelo que mejor aproxima el comportamiento de los geles es el modelo de Steinwachs, seguido por el de Holzapfel, siendo el Neo-Hookeano el que peores resultados ofrece. Por último, se ha desarrollado una sub-rutina UMAT que implementa el modelo de Steinwachs en el entorno de Análisis por Elementos Finitos Simulia ABAQUS, con el objetivo de ser empleada en el futuro en la resolución de problemas de mayor relevancia en el contexto de TFM. Palabras clave: Mecanobiología,Microscopía de Fuerza de Tracción,Hiperelasticidad ,Mecánica no-lineal,Método de Elementos Finitos,MATLAB,ABAQUS,Trabajo de Fin de Máster. Abstract This master’s thesis is devoted to the study of the viability of the application of three hyperelastic constitutive models for the characterization of collagen hydrogels. The considered models are: Neo-Hookean model, Holzapfel model and Steinwachs model. The gels to be characterized are used in mechanobiology to simulate the characteristics of the extracellular matrix, being of great importance in the framework of Traction Force Microscopy (TFM). In such a study, the implementation of the models has been carried out in MATLAB, and using the tools offered by the software, it has been performed a least squares fitting of the experimental curves, sent by collaborators from the KU Leuven University, to the mentioned models. The model that best approximates the behaviour of the gels is the Steinwachs model, followed by Holzapfel’s, with the Neo-Hookean model offering the worst results. Finally, a UMAT subroutine aimed for the implementation of the Steinwachs model in the Finite Element Analysis environment Simulia ABAQUS, has been developed, which could be used in the future in relevant problems in the framework of TFM. Keywords: Mechanobiology,Traction Force Microscopy,Hiperelasticity ,Non-linear mechanics, Finite Element Method,MATLAB,ABAQUS,Master’s Thesis. Agradecimientos Querría agradecer la supervisión de este trabajo al Profesor José Antonio Sanz Herrera, tutor del mismo, cuyos conocimientos y opiniones han resultado indispensables para guiar mi camino en el proceso. Asimismo, querría agradecer al equipo de investigación en Mecanobiología e Ingeniería de Tejidos dirigido por el Profesor Hans Van Oosterwyck en la Universidad KU Leuven los datos experimentales en torno a los cuales han girado la mayoría de actividades realizadas en este trabajo. movimiento de los beads inmersos en el substrato en lo que se conoce como 3D TFM, y constituye un procedimiento relativamente reciente [4], habiéndose practicado con menor frecuencia que la modalidad bidimensional. Esto se debe, no exclusivamente al hecho de que el problema implicado se manifieste de mayor complejidad que en 2D TFM, sino además a que esta última tiene fundamento en la observación de que a menudo las células se vuelven “planas” en adhesiones maduras dentro de substratos elásticos. No obstante, muchos tipos celulares requieren la reconstrucción de las complejas características biofísicas y bioquímicas de una matriz extracelular tridimensional para poder mostrar un comportamiento fisiológicamente realista en el interior de ambientes in vitro [5]. Asimismo, cabe discernir entre procedimientos que contemplan no-linealidades de origen geométrico (deformaciones finitas, ε∼O(1)) y aquellos que no (deformaciones infinitesimales, ε1). Dentro del régimen lineal, la aproximación de deformaciones infinitesimales brinda la posibilidad de usar funciones de Green en la resolución del problema mecánico, dando lugar a rápidos y eficientes algoritmos para la reconstrucción de las tracciones. Sin embargo, avances recientes en TFM muestran la capacidad de las células de ejercer deformaciones con una magnitud de hasta el 40 %, requiriéndose un enfoque basado en la teoría de grandes deformaciones que contemple tales niveles de deformación [6]. A este respecto, la resolución del problema mecánico habrá de enmarcarse en el contexto de una discretización por elementos finitos, no pudiéndose utilizar las funciones de Green para ello. Ambos procedimientos se hallan inscritos en el marco del método inverso de recuperación de tracciones. Esencialmente, este método se basa en la obtención de las tracciones a través de la resolución de un problema de minimización que implique el equilibrio de fuerzas internas y de volumen (∇ · σ+f= 0), así como determinadas restricciones, empleando como dato el campo de desplazamientos medido experimentalmente. De esta forma, el campo de tracciones resultante no se corresponderá directamente con el campo de desplazamientos medido. Existe, no obstante, una alternativa a la resolución del problema inverso: el método directo; el cual se basa en la obtención del tensor de deformaciones correspondiente al campo de desplazamientos medido (derivándolo) para, mediante la ley constitutiva del material, obtener el tensor de tensiones asociado. Ambos métodos se describen con más detenimiento en las secciones 1.2.1 y 1.2.2. La figura 1.1 ilustra una configuración típica en 3D TFM. La superficie azul corresponde a una célula de tumor de mama inmersa en una matriz de colágeno de tipo I, y las barras coloreadas representan los desplazamientos correspondientes a los beads causados por la relajación de la tracción celular [7]. Figura 1.1: Una célula de tumor de mama (azul) inmersa en una matriz de colágeno tipo I. 2 1.2.1. Método directo El método directo consiste fundamentalmente en la reconstrucción directa de las tracciones celulares a partir del campo de desplazamientos medido en laboratorio. La diferenciación de este último permite obtener el tensor gradiente de deformación correspondiente, que es entonces aplicado a la ley constitutiva del material para obtener el tensor de tensiones de Cauchy. Las tracciones buscadas se relacionan con dicho tensor a través del vector normal exterior genérico de la superficie envolvente de la célula. Un flujo de trabajo típico para este método podría ser: 1. Obtención del campo de desplazamientos discreto, esto es, los desplazamientos de cada bead. 2. Interpolación mediante funciones de forma del campo de desplazamientos en el dominio considerado. Esto requerirá la discretización del mismo en elementos finitos, pero aliviará bastante la generación de ruido a la hora de practicar las derivadas sobre el campo de desplazamientos, la cual puede devenir enorme en la aproximación numérica de derivadas. 3. Diferenciación del campo de desplazamientos y obtención del tensor gradiente de deformación compatible a través de la expresión F=1+ [∇⊗u]T. 4. Obtención del tensor de tensiones de Cauchy, σ, a través de la relación constitutiva del substrato, que ha de depender en última instancia de F. 5. Recuperación de las tracciones en la interfaz célula-substrato, la cual se relaciona con el tensor de tensiones de Cauchy a través del vector normal a la superficie de la célula, t=σ·n. Debido al carácter “ruidoso” de los desplazamientos medidos, el método directo presenta potenciales inestabilidades en su aplicación, esto es, gran sensibilidad de las fuerzas reconstruidas ante pequeños cambios en los desplazamientos medidos. Por tanto, el procedimiento requiere de la disposición de técnicas sofisticadas de procesamiento de imágenes, así como de la captura de imágenes 3D de gran resolución. 1.2.2. Método inverso Este enfoque busca una solución regularizada matemáticamente consistente del campo de tracciones, que en última instancia implique a las relaciones constitutivas y cinemáticas (de compatibilidad) entre desplazamientos y deformaciones en una región grande alrededor de la célula. El objetivo del procedimiento es pues el de obtener una solución al campo de tracciones a través de un problema de minimización con restricciones que contemple como parte de la función objetivo la ecuación de equilibrio que relaciona las tracciones con el campo de desplazamientos medido. En el caso de deformaciones infinitesimales, dicha relación se expresa a través de la integral de convolución de la función de Green correspondiente con las tracciones buscadas, que en la práctica es discretizada y por tanto pasa a ser una operación matricial, u(x) = ZG(X,x)t(x)dx−→ u=G∗t(1.1) donde ues el campo de desplazamientos medido, G∗tes el campo de desplazamientos predicho por el campo de tracciones estimado (Gla función de Green, tel campo de tracciones buscado), es decir, la integral de convolución discretizada; y xyXlas coordenadas de cualquier punto del substrato en configuración deformada e indeformada, respectivamente. Debido a la naturaleza tensorial de la teoría de la elasticidad, la función de Green Gconstituye una matriz 3×3la cual contempla desplazamientos del substrato en tres dimensiones. De este modo, la solución de (1.1) se obtiene de resolver un problema de minimización por mínimos cuadrados con regularización de la función objetivo (u−Gt). La regularización constituye un enfoque heurístico que tiene por objetivo limitar la magnitud de las tracciones resultantes de resolver dicho problema, que puede alcanzar valores muy grandes debido al ruido presente en las mediciones del campo de desplazamientos. Atendiendo a tal restricción, se plantea el problema de minimización con regularización de Tikhonov de orden cero como [3], 3 t= arg m´ın th[Gt −u]2+λ2t2i(1.2) en donde la elección del parámetro de regularización λpuede basarse en la teoría Bayesiana [8] o realizarse empíricamente [9]. Por tanto, la diferencia fundamental del método inverso con respecto al directo es que el campo de tracciones resultante no guarda relación directa con el campo de desplazamientos medido. En caso de que las deformaciones sean grandes, el problema de optimización descrito deberá plantearse en el contexto de una discretización por elementos finitos u otros métodos afines. En todo caso, el método inverso representa pues un problema de gran complejidad computacional en relación al método directo, pero proporciona resultados significativamente más fidedignos que este último. Por ello, para conjuntos de datos de calidad estándar, resulta más conveniente aplicar el método inverso. 1.3. Objetivo El objetivo del presente trabajo consiste en verificar la viabilidad de los modelos hiperelásticos Neo-Hookeano, de Holzapfel y de Steinwachs en la caracterización de las propiedades mecánicas de unos hidrogeles de colágeno tipo I. Para ello, se pretende implementar dichas leyes constitutivas en MATLAB y hacer uso de las herramientas de optimización que ofrece tal programa a la hora de ajustar los modelos a datos experimentales del comportamiento de los geles enviados por colaboradores en la Universidad KU Leuven. Por último, se quiere desarrollar una sub-rutina de usuario que permita implementar el modelo material de Steinwachs en el programa de análisis por elementos finitos Simulia ABAQUS. 4 Capítulo 2 Hiperelasticidad En este capítulo se desarrolla el núcleo teórico del trabajo. Se comienza realizando una breve descripción del contexto teórico en el que se enmarcan las subsiguientes secciones, esto es, el compuesto por la formulación de las magnitudes esenciales de tensión y deformación en la elasticidad no lineal, y la teoría de leyes constitutivas hiperelásticas. A continuación, se describen de forma detallada las características de las tres leyes hiperelásticas consideradas en este trabajo. 2.1. Marco teórico Aunque el análisis lineal de problemas de elasticidad sea, en general, de gran utilidad -además de relativamente simple-, su alcance es limitado. Las condiciones requeridas para que su aplicación sea viable no son satisfechas en muchas de las circunstancias dadas en la práctica, como aquellas propias de la mecanobiología. En efecto, cualquier situación que involucre la respuesta en deformación de una matriz extracelular ante un determinado estado de cargas -como, por ejemplo, un gel de colágenono podrá ser tratada desde un enfoque lineal. Esto es cierto para virtualmente todo tejido corporal o elemento mecánico biológico. Existen, pues, dos clases esenciales de no-linealidades: las no-linealidades de origen material, y aquellas de origen geométrico. La primera se debe a que la propia ley constitutiva del material no se comporta de forma lineal, es decir, la relación tensión-deformación no es una línea recta, al contrario que aquella de un material elástico lineal (σ=C:ε). La segunda, en cambio, se origina cuando los desplazamientos producidos en el proceso de deformación son lo suficientemente grandes como para que las configuraciones deformada e indeformada del cuerpo en cuestión no puedan ser consideradas indistintas. Todas las definiciones y formulaciones que se hallan en la presente sección han sido extraídas fundamentalmente de los libros de Holzapfel [10] y Bonet yWood [11]. 2.1.1. Breve introducción a las medidas de tensión-deformación en elasticidad no lineal Contrariamente a como se hace en la teoría de deformaciones infinitesimales, en la teoría de grandes deformaciones (elasticidad geométricamente no lineal) se considera la distinción entre las configuraciones deformada e indeformada del cuerpo en deformación, debido a que los desplazamientos que las relacionan son, como mínimo, del orden de magnitud de las dimensiones originales del problema (configuración indeformada). De este modo, la posición instantánea de un punto cualquiera del medio puede expresarse a partir de las coordenadas iniciales como, x(X, t) = X+u(X, t)→u(X, t) = O(X)(2.1) 5 donde ues el vector de desplazamientos. Esta diferenciación entre xyXda lugar a dos tipos de formulaciones distintas de las magnitudes y ecuaciones de la mecánica no lineal: la formulación material o Lagrangiana, en la cual se realiza el seguimiento de la evolución de las propiedades de una partícula que en el instante inicial se hallaba en un punto determinado de la configuración indeformada; y la formulación espacial o Euleriana, en la que la atención se centra más bien en la evolución de las magnitudes en cualquier punto del dominio definido, siendo estos puntos ocupados por diferentes partículas en diferentes instantes. Así pues, la primera se enfoca en la partícula, mientras que la segunda lo hace en el dominio en el que esa partícula se halla. Siendo θ una magnitud cualquiera, ambas configuraciones se “encuentran” a través de la siguiente relación: θ(X, t) = θ(x(X, t), t)(2.2) Se define el tensor gradiente de deformación Fcomo la magnitud que permite expresar la posición relativa espacial entre dos partículas vecinas después de la deformación en términos de su posición relativa material previa a la deformación, F=∂x ∂XFiJ =∂xi XJ∀i, J = 1,2,3(2.3) En este sentido, Fpermite transformar vectores expresados en la configuración material, o de referencia, en vectores expresados en la configuración espacial, o actual; y por ello es conocido como un tensor bi-punto. Esto quiere decir que no se halla expresado en ninguna de las dos configuraciones, sino “en medio” de ambas. De ahí la diferenciación entre subíndices en minúsculas (espacial) y mayúsculas (material). Del mismo modo pueden definirse otras magnitudes de deformación que operen exclusivamente sobre elementos materiales o espaciales: Materiales Tensor de Cauchy por la derecha,C=FTF Tensor de deformaciones de Green,E= [C−1]/2 Espaciales Tensor de Cauchy por la izquierda,b=F F T Tensor de deformaciones de Almansi,e= [1−b−1]/2 En cuanto a las magnitudes de tensión, aunque éstas se describen con mayor profundidad junto con la introducción de la noción de función de densidad de energía de deformación en la siguiente sección, puede anticiparse cuáles son de naturaleza material, bi-punto y espacial: Material.Segundo tensor de Piola-Kirchhoff (S). Bi-punto.Primer tensor de Piola-Kirchhoff (P). Espacial.Tensor de Cauchy (σ) y Tensor de Kirchhoff (τ). Resulta conveniente definir un operador a través del cual puedan expresarse las relaciones entre las distintas medidas de tensión descritas por los tensores presentados previamente. Se define pues el operador de empuje opush-forward como aquel que permite transformar magnitudes materiales en espaciales, φ∗[]; y el de tirón opull-back como aquel que permite transformar magnitudes espaciales en materiales, φ∗[]. Ambos son mutuamente inversos, y la operación implicada por ellos dependerá de la magnitud sobre la cual se practique. De este modo, las relaciones entre las distintas medidas de tensión a través del operador de empuje son las siguientes: τ=Jσ=φ∗[P] = P F T(2.4) τ=Jσ=φ∗[S] = F SF T(2.5) 6 2.1.2. Hiperleasticidad Se define la elasticidad como la cualidad de aquel material cuya relación constitutiva es únicamente dependiente del estado actual de deformación. En tal sentido, la hiperelasticidad es definida como la cualidad propia de aquellos materiales elásticos cuya relación constitutiva deriva directamente de una función escalar que relaciona la energía de deformación con el tensor gradiente de deformación, denominada función densidad de energía de deformación opotencial elástico. Esto es posible debido a que el trabajo realizado por las tensiones en el proceso de deformación depende exclusivamente de los instantes inicial y final, no por tanto del camino intermedio. Esta función se define como, Ψ(F(X),X) = Zt to P(F(X),X):˙ Fdt −→ ˙ Ψ = P:˙ F(2.6) donde la notación ˙ [] indica derivada con respecto al tiempo. Toda función densidad de energía debe ser objetiva, esto es, debe permanecer invariante ante movimientos de sólido rígido. Teniendo en cuenta que el tensor gradiente de deformación puede ser descompuesto en la combinación de una rotación y un estiramiento (descomposición polar), F=RU, el potencial elástico ha de depender únicamente de la componente de estiramiento U. Por una cuestión de conveniencia, el potencial se suele expresar en función de U2=C, permitiendo así la reinterpretación de la expresión dada en (2.6) como, Ψ(C(X),X) = Zt to S(C(X),X):˙ Edt −→ ˙ Ψ = S:˙ E(2.7) Los pares Py˙ FySy˙ Eson conjuados energéticos, y permiten definir la relación constitutiva en las configuraciones mixta y material, respectivamente, de la siguiente forma, P(F(X),X) = ∂Ψ(F(X),X) ∂FPiJ =∂Ψ ∂FiJ (2.8) S(C(X),X) = ∂Ψ(C(X),X) ∂ESIJ =∂Ψ ∂EIJ (2.9) Asimismo, en el marco de un material isótropo, el potencial elástico debe ser independiente de cualquier dirección en la configuración material, i. e., ha de depender únicamente de los tres primeros invariantes del tensor C, Ψ(C(X),X) = Ψ(I1, I2, I3,X)(2.10) en donde los invariantes de Cson definidos como, I1=tr[C] = C: 1 (2.11) I2=tr[C2] = C:C(2.12) I3= det[C] = det[F]2=J2(2.13) Las relaciones dadas en (2.8) y (2.9) serán, en general, no lineales. En el contexto de un esquema numérico iterativo -por ejemplo, un esquema de Newton-Raphson- aparecerán linealizadas con respecto a un incremento uen la configuración correspondiente, presentándose así dichas relaciones en la siguiente forma, DP[u] = A:DF[u]DS[u] = C:DE[u](2.14) 7 en donde Df[u]denota derivada direccional de fen la dirección de u.AyCson tensores de cuarto orden que reciben el nombre de Primer ySegundo Tensor de Elasticidad, respectivamente, y vienen descritos por, A=∂P ∂FAiJkL =∂PiJ ∂FkL (2.15) C=∂S ∂ECIJKL =∂SIJ ∂EKL (2.16) Ambos tensores, debido a que son consecuencia de derivar dos veces una función (potencial elástico) con respecto a la misma cantidad (tensor de deformación), poseen simetría mayor, AiJkL =AkLiJ CIJKL =CKLIJ (2.17) antes bien, sólo el segundo tensor de elasticidad presenta simetrías menores, ya que es consecuencia de derivar una cantidad simétrica (segundo tensor de Piola-Kirchoff) con respecto a otra magnitud simétrica (tensor de deformación de Green-Lagrange), CIJKL =CJIKL CIJKL =CIJLK (2.18) 2.2. Leyes hiperelásticas estudiadas A continuación se describen en detalle las tres leyes hiperelásticas consideradas en este trabajo: la ley Neo-Hookeana, la ley de Holzapfel y la ley de Steinwachs. 2.2.1. Neo-Hookeana El modelo Neo-Hookeano constituye una extensión del comportamiento descrito por la ley de Hooke para el caso de grandes deformaciones [12]. A diferencia con el sólido Hookeano, cuya ley de comportamiento es enteramente lineal, la relación tensión-deformación del sólido Neo-Hookeano es virtualmente lineal antes de cruzar el umbral de pequeñas deformaciones, más allá del cual el comportamiento deviene no lineal. La tensión alcanza eventualmente un plateau debido a la liberación de energía de deformación en forma de calor, para luego volver a rigidizarse. Fue concebido en origen para describir el comportamiento de materiales plásticos y similares a la goma, de carácter incompresible, es decir, cuyo volumen permanece invariable durante el proceso de deformación (J= 1), y por tanto su función densidad de energía de deformación no depende del tercer invariante. No obstante, en ocasiones conviene llevar a cabo la concepción de modelos compresibles, ya sea porque la implementación numérica de tal modelo presente dificultades [13], o porque el propio material de estudio sea de naturaleza compresible, como es el caso del considerado en este trabajo, un hidrogel de colágeno de tipo I cuasi-incompresible [14, 15]. A este respecto, la función potencial elástico del material ahora habrá de depender del tercer invariante , o lo que es igual, del determinante de F(J). Existen diversas formulaciones alternativas del potencial elástico para materiales Neo-Hookeanos compresibles, habiéndose elegido para este trabajo la siguiente [13], Ψ(I1, J) = µ 2[I1−1 : 1 −2 log J] + κ 2[J−1]2(2.19) siendo µyκlos parámetros de Lamé, con unidades de tensión. Se verifica que cuando no existe deformación, F=1, el potencial elástico adquiere valor cero. La ley constitutiva en configuración mixta viene descrita por, P=∂Ψ(I1, J) ∂F=∂Ψ ∂I1 ∂I1 ∂F+∂Ψ ∂J ∂J ∂F(2.20) 8 siendo, ∂Ψ ∂I1 =µ 2 ∂I1 ∂F= 2F(2.21) ∂Ψ ∂J =−µ1 J+κ[J−1] ∂J ∂F=J·[F−1]T=cof[F](2.22) denotando cof[F]la matriz de cofactores de F. 2.2.2. Holzapfel En 1998, Holzapfel yWeizsäcker [16] propusieron un modelo desacoplado del potencial elástico para arterias de tipo elástico ymuscular, Ψ=Ψiso + Ψaniso (2.23) El término Ψiso representa la contribución isótropa de la sustancia base no colágena, probablemente por la acción de la elastina en la pared vascular, y que puede corresponder a una función densidad de energía de deformación isótropa cualquiera. En este sentido, Ψiso será una función de los tres primeros invariantes, y en el contexto de este trabajo se ha elegido como la descrita en la sección 2.2.1. Por otro lado, Ψaniso representa la contribución anisótropa de la familia de fibras colágenas presentes en la pared vascular en una dirección determinada. Ψaniso será pues una función del invariante cuarto, el cual se define como, I4=MTCM (2.24) en donde Mcorresponde al vector unitario que apunta en la dirección paralela a las fibras. La manera en la que Ψaniso depende de I4es a través de una función exponencial propuesta originalmente por Fung et al. [17], Ψaniso(I4) = k1 k2hexp hk2[I4−1]2i−1i(2.25) y que contempla la fuerte rigidización del tejido a cargas más elevadas debida casi enteramente a la acción de las fibras colágenas. k1es un parámetro material con unidades de tensión y k2es un parámetro adimensional. Si se consideran más de una familia de fibras, la contribución total se obtendrá como el promedio de las contribuciones de toda ellas, pasando entonces (2.25) a escribirse como, Ψaniso(I4) = k1 k2 1 n n X i=1 hexp hk2[I4i−1]2i−1i(2.26) siendo nel número total de familias de fibras (o direcciones) existentes. Se verifica que para F=1, esto es, en ausencia de deformaciones, I4= 1, y Ψaniso deviene nulo. No obstante, el gel de colágeno considerado en este trabajo se comporta de forma isótropa. Esto implica que las fibras embebidas en el mismo han de estar homogéneamente distribuidas en todas direcciones, esto es, isótropamente distribuidas. Debido a esto, para obtener la ley constitutiva correspondiente, el término Ψaniso debe constituir el promedio de la contribución de todas las fibras existentes en todas direcciones. Así pues, la ley constitutiva en configuración mixta queda expresada de la siguiente forma, P=∂Ψ(I1, J, I4) ∂FΩ =∂Ψiso(I1, J) ∂F+∂Ψaniso(I4) ∂FΩ = =∂Ψiso ∂I1 ∂I1 ∂F+∂Ψiso ∂J ∂J ∂F+∂Ψaniso ∂I4 ∂I4 ∂FΩ (2.27) 9 donde los términos correspondientes a la parte isótropa se obtienen de forma idéntica a lo indicado en (2.21) y (2.22), mientras que aquellos correspondientes a la parte anisótropa se obtienen como, ∂Ψaniso ∂I4 = 2k1[I4−1] exp k2[I4−1]2∂I4 ∂F= 2[F M]⊗M(2.28) denotando a⊗bel producto de Kronecker entre ambos vectores. La notación h·iΩindica el promedio de la magnitud entre paréntesis angulares considerándose toda dirección contenida en el ángulo sólido Ω= [θ, φ], donde θyφson las variables independientes angulares en coordenadas esféricas; y Mun vector unitario cualquiera definido como M= [sin θcos φ, sin θsin φ, cos θ]T. Se ha elegido realizar el promedio integrando la magnitud en cuestión alrededor de la esfera unitaria, para ello empleando la cuadratura de Gauss-Legendre [18], ∂Ψaniso ∂I4 ∂I4 ∂FΩ =1 4πZ2π 0Zπ 0 ∂Ψaniso ∂I4 ∂I4 ∂F(θ, φ) sin θdθdφ ≈ ≈1 4π π m 2m X j=1 m X i=1 wi∂Ψaniso ∂I4 ∂I4 ∂F(θi, φj)(2.29) siendo mel orden de integración. El conjunto {θi}es escogido de modo que {cos (θi)}y{wi}sean, respectivamente, los nodos y pesos de Gauss-Legendre en el intervalo [−1,1]. Los puntos φjse distribuyen uniformemente en el intervalo [0,2π]a través de la siguiente expresión, φj=j−1 2π m(2.30) 2.2.3. Steinwachs Steinwachs et al. [19] proponen un modelo para redes de biopolímeros -como el colágeno-, basado en la combinación de una descripción micromecánica, es decir, centrada en el comportamiento individual de una fibra en particular; y de una descripción de medio continuo. Teniendo en cuenta una escala espacial suficientemente pequeña correspondiente a un segmento de fibra, la deformación local de dicho segmento es considerada no-afín, esto es, independiente de la deformación del resto del volumen. Sin embargo, para escalas mayores a la distancia de interconexión típica entre fibras, la deformación de la fibra se aproxima a la deformación macroscópica λ, dependiendo de la orientación de la fibra y la deformación practicada. De este modo, se considera que las deformaciones pasan a ser afines para un volumen de material suficientemente grande. Análogamente a lo descrito en la sección 2.2.2, la relación constitutiva macroscópica para este modelo se obtiene promediando en todas direcciones aquella propia de una fibra. Se considera que una fibra puede encontrarse en tres estados distintos de deformación: compresión, ante la cual la fibra se vuelve inestable (pandea) y presenta entonces una rigidez que decae exponencialmente hasta cero; enderezamiento, proceso intermedio en el que la fibra recupera su longitud natural después de haber pandeado y ante el cual presenta una rigidez constante; y estiramiento, proceso en el que la fibra se deforma axialmente más allá de su longitud natural y experimenta una rigidización exponencial. El potencial elástico se define a través de la rigidez, es decir, de su segunda derivada con respecto a la deformación macroscópica λ, Ψ00(λ) = κo   exp [λ/do]∀λ < 0 1∀0≤λ<λs exp [[λ−λs]/ds]∀λ≥λs λ=kF eΩk−1(2.31) en donde κoes un parámetro material con dimensiones de tensión, does el parámetro adimensional de pandeo, dses el parámetro adimensional de estiramiento y λses el parámetro adimensional que controla la extensión del intervalo de comportamiento lineal, i. e., de rigidez constante. La deformación λde una fibra se define como la variación de longitud de la misma con respecto a su 10 longitud natural, y eΩes el vector unitario que describe la orientación de la fibra. Se comprueba que el potencial elástico presenta un valor nulo cuando F=1, puesto que en ese caso λ= 0 (no existe deformación alguna). La tensión mecánica de la fibra es definida por la siguiente expresión, Ψ0(λ) = Zλ 0 Ψ00(λ)dλ =κo   do·[exp [λ/do]−1] ∀λ < 0 λ∀0≤λ<λs λs+ds·[exp [[λ−λs]/ds]−1] ∀λ≥λs (2.32) indicando la integración desde cero que el material no se encuentra pretensado. La ecuación constitutiva en configuración mixta se escribe, P=∂Ψ(λ(F)) ∂FΩ =∂Ψ ∂λ ∂λ ∂FΩ =Ψ0(λ)[F eΩ]⊗eΩ kF eΩkΩ (2.33) La Figura 2.1 muestra de forma ilustrativa el comportamiento de los tres modelos en el caso de una extensión uniaxial (movimiento impedido en las direcciones yyz), caracterizada por el siguiente gradiente de deformación, F=  α0 0 0 1 0 0 0 1  (2.34) Los parámetros de los modelos se han elegido para que las curvas resultantes sean similares entre sí. Para la curva Neo-Hookeana, el conjunto de valores es: (µ, κ) = (1.5,20) [Pa]; para la curva de Holzapfel, el conjunto de valores es: (µ, κ, k1, k2) = (1.5Pa,20 Pa,0.85 Pa,0.06); y para la curva de Steinwachs el conjunto de valores es: (κo, ds, λs, do) = (200 Pa,0.5,0.5,0.005). Cabe destacar que, de las tres curvas, aquella producida por el modelo de Steinwachs es la única que presenta rigidez nula en la región de compresión. Esto es, naturalmente, debido a la propia definición del potencial elástico del modelo, la cual establece que las fibras pandean rápidamente después de ser sometidas a compresión. El modelo de Holzapfel mostraría un comportamiento similar en caso de no contemplar la aportación de la matriz, esto es, la contribución isótropa del potencial Neo- Hookeano. Este fenómeno, tal y como se verá en secciones posteriores, es susceptible de entrañar numerosos problemas de carácter numérico. 0 0.5 1 1.5 2 2.5 3 α -200 -100 0 100 200 300 400 σ11 [Pa] Neo-Hook Holzapfel Steinwachs Figura 2.1: Tensión de Cauchy frente a deformación en extensión uniaxial. 11 0 0.5 1 1.5 2 2.5 3 α -50 0 50 100 150 200 250 300 σ11 [Pa] do = 0.005 do = 0.1 do = 0.5 Figura 3.13: Tensión de Cauchy en extensión uniaxial para diferentes valores de do(modelo de Steinwachs). -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 γ -150 -100 -50 0 50 100 150 σ12 [Pa] do = 0.005 do = 0.1 do = 0.5 Figura 3.14: Tensión de Cauchy en cizalla simple para diferentes valores de do(modelo de Steinwachs). 3.3. Procedimiento de ajuste de los datos experimentales El ajuste de los parámetros de los modelos a datos experimentales se ha llevado a cabo a través del empleo de la función de MATLAB lsqnonlin [25], la cual resuelve problemas no-lineales de ajuste de curvas por mínimos cuadrados de la forma, m´ın xkf(x)k2 2= m´ın x(f1(x)2+f2(x)2+···+fn(x)2)(3.8) El procedimiento seguido se resume en los siguientes pasos: 1. Se ajustan heurísticamente los parámetros para obtener una curva relativamente similar a la experimental. Los valores encontrados se almacenan en un vector param0. Estos valores de partida favorecerán la convergencia del algoritmo, evitando en la medida de lo posible la aparición de mínimos locales. 2. Se selecciona el rango de abscsisas en el que se quiere ajustar. Entonces se almacenan las abscsisas y ordenadas experimentales en dos vectores xdata yydata, respectivamente. Las abscisas corresponden a los valores de deformación del experimento, así como las ordenadas a los valores de tensión correspondientes. 3. Se define la función correspondiente al modelo hiperelástico modelo(param,xdata), de modo que acepte como variables de entrada el vector de parámetros param y el vector de las abscisas experimentales xdata, devolviendo las ordenadas predichas por la ley constitutiva en cuestión. 4. Se define una función anónima que compute la diferencia entre las ordenadas predichas por el modelo y las experimentales: fun = @(param)modelo(param,xdata) - ydata. 5. Opcionalmente se definen dos vectores lb yub que contienen, respectivamente, cotas inferiores y superiores de los valores que pueden alcanzar los parámetros. 6. Se arranca el cálculo: param = lsqnonlin(fun,param0,lb,ub,options). La variable de entrada options es opcional, y permite la especificación de algoritmos de resolución no estándar, así como de otras opciones. En este trabajo se ha omitido, utilizándose por tanto el conjunto de opciones definidas por defecto. 18 3.3.1. Obtención de las curvas de contracción vertical El ajuste de curvas de contracción vertical entraña el problema de que el conjunto de ordenadas se corresponde con una componente del tensor gradiente de deformación, siendo entonces este último desconocido y no pudiendo por tanto ser utilizado como variable de entrada directa en la ley constitutiva. En efecto, el tensor gradiente de deformación que caracteriza un estado de contracción vertical ante extensión horizontal se escribe como, F=  λh0 0 0 1 0 0 0 λv (3.9) donde λhes la extensión horizontal aplicada (abscisas) y λves la contracción vertical resultante de aplicar dicha extensión (ordenadas). En este sentido, λvse obtiene de aplicar la condición de libre movimiento en la dirección vertical, esto es, σ33 = 0. A continuación se procede a explicar la obtención de la contracción vertical λvpara cada modelo. Para aproximar las integrales de superficie se ha utilizado el método de cuadratura de Gauss-Legendre descrito en la sección 2.2.2. Neo-Hookeano El tensor de tensiones de Cauchy para el modelo Neo-Hookeano recibe la siguiente expresión, σ=1 JµF F T+1 J−µ1 J+κ[J−1]·cof[F]FT(3.10) siendo, J=λhλv,cof[F] =   λv0 0 0λhλv0 0 0 λh (3.11) A través de la aplicación de la condición σ33 = 0 se tiene, [µ+κλ2 h]λ2 v−κλhλv−µ= 0 (3.12) Reteniendo la única solución compatible (positiva), se obtiene la siguiente expresión explícita para λv, λv= κλh+qκ2λ2 h+ 4µ[µ+κλ2 h] 2 [µ+κλ2 h](3.13) Holzapfel Debido al carácter fuertemente no-lineal del modelo de Holzapfel -y de Steinwachs-, no es posible obtener una expresión explícita para la contracción vertical. Por tanto, resulta necesario el planteamiento de un esquema iterativo para aproximar la solución, habiéndose elegido en este caso el método de Newton-Raphson-Kantorovich. Este último se enuncia de la siguiente manera, f(x)=0 −→ xk+1 =xk−f(xk) f0(xk)(3.14) particularizándose en este caso como, σ33(λv)=0 −→ λk+1 v=λk v−σ33(λk v) σ0 33(λk v)(3.15) donde el índice kdenota la iteración actual. Es preciso pues, para cada k+ 1, partir de un valor inicial (semilla) en kque arranque el proceso iterativo. En este caso, para incrementos 19 suficientemente pequeños en λh(|λn+1 h−λn h|  1), se asume que λk v|λn+1 h≈λv|λn hal principio de la iteración, teniéndose en cuenta que para λh= 1,λv= 1. Nótese que el índice nindica la estación actual, esto es, la posición actual en la variable independiente λh, y por tanto es distinto de k. En este sentido, una estación ncontiene múltiples iteraciones k. El tensor de tensiones de Cauchy se obtiene mediante la siguiente expresión, σ=1 JµF F T+1 J−µ1 J+κ[J−1]·cof[F]FT+ +1 J1 4πZ2π 0Zπ 0 2k1[I4−1] exp k2[I4−1]2·2[F M]⊗M·sin θdθdφ·FT(3.16) Imponiendo σ33 = 0, σ33(λv) = 1 λh µλv−µ1 J+κ[J−1]+ +1 λh 1 4πZ2π 0Zπ 0 2k1[I4−1] exp [k2[I4−1]2]·2λvM2 3·sin θdθdφ = =1 λhf1+f2+1 4πZ2π 0Zπ 0 2f3M2 3·sin θdθdφ≈ ≈1 λh f1+f2+1 4π π m 2m X j=1 m X i=1 wi2f3M2 3(θi, φj) = 0 (3.17) donde, f1=µλv, f2=−µ1 J+κ[J−1]λh, f3= 2k1[I4−1] exp k2[I4−1]2λv(3.18) La derivada de σ33 con respecto a λv, σ0 33(λv) = 1 λhf0 1+f0 2+1 4πZ2π 0Zπ 0 2f0 3M2 3·sin θdθdφ≈ ≈1 λh f0 1+f0 2+1 4π π m 2m X j=1 m X i=1 wi2f0 3M2 3(θi, φj) (3.19) donde, f0 1=µ, f0 2=µ1 λ2 v +κλ2 h, f0 3=22k1+ 4k1k2[I4−1]2M2 3λ2 v+ 2k1[I4−1]·exp k2[I4−1]2(3.20) Steinwachs Análogamente al caso del modelo de Holzapfel, la obtención de la contracción vertical para el modelo de Steinwachs se realiza mediante el planteamiento de un esquema iterativo basado en el método de Newton-Raphson-Kantorovich, descrito por (3.15). El tensor de tensiones de Cauchy se obtiene como, 20 σ=1 J1 4πZ2π 0Zπ 0 Ψ0(λ)[F eΩ]⊗eΩ kF eΩksin θdθdφ·FT(3.21) Imponiendo σ33 = 0, σ33(λv) = 1 λh1 4πZ2π 0Zπ 0 Ψ0(λ)1 kF eΩkλvcos2θsin θdθdφ= =1 λh1 4πZ2π 0Zπ 0 f·cos2θsin θdθdφ≈ ≈1 λh  1 4π π m 2m X j=1 m X i=1 wif(θi, φj)·cos2θi = 0 (3.22) donde, f= Ψ0(λ)1 kF eΩkλv(3.23) La derivada de σ33 con respecto a λv, σ0 33(λv) = 1 λh1 4πZ2π 0Zπ 0 f0·cos2θsin θdθdφ≈ ≈1 λh  1 4π π m 2m X j=1 m X i=1 wif0(θi, φj)·cos2θi (3.24) donde, f0= Ψ00(λ)dλ dλv 1 kF eΩkλv+ Ψ0(λ)"kF eΩk− dkF eΩk dλvλv kF eΩk2#(3.25) dλ dλv =dkF eΩk dλv =λvcos2θ kF eΩk(3.26) 21 Capítulo 4 Ajuste experimental de los parámetros gobernantes de las leyes En este capítulo se presentan los principales resultados obtenidos de ajustar los parámetros gobernantes de las leyes constitutivas descritas en la sección 2.2 a las curvas experimentales recibidas mediante el procedimiento descrito en la sección 3.3. Estas curvas corresponden a dos experimentos en concreto: cizalla simple y contracción vertical ante extensión horizontal (efecto Poisson); practicados ambos sobre dos hidrogeles de colágeno con sendas concentraciones de fibras: 0.6 y 1.2 mg/ml. La caracterización reológica de los hidrogeles (cizalla simple) ha sido llevada a cabo por Mar Cóndor Salgado, integrante del equipo de investigación en Mecanobiología e Ingeniería de Tejidos dirigido por el Profesor Hans Van Oosterwyck en la Universidad KU Leuven, siguiendo las pruebas de Steinwachs et al. [19]. Los resultados correspondientes a la contracción vertical han sido extrapolados de los datos obtenidos por Ben Fabry, coautor del artículo citado en líneas anteriores. 4.1. Datos experimentales recibidos 4.1.1. Cizalla simple El gel se somete a la acción de un reómetro rotacional de cono y plato con el fin de medir la relación tensión-deformación en un estado de deformación de cizalla simple. Para ello, el gel se deposita en un plato horizontal y a su vez se sitúa un cono poco profundo encima del mismo. La polimerización del gel se efectúa en el interior de la configuración. Entonces, el reómetro aplica un estado de deformación de modo que el potencial elástico del gel sólo depende de la deformación de cizalla. Básicamente, el plato es rotado y el torque transmitido al cono es medido, a partir del cual puede obtenerse la tensión de cizalla σ12. El tensor gradiente de deformación que caracteriza el experimento recibe la expresión, F=  1γ0 0 1 0 0 0 1  (4.1) donde γrepresenta la deformación de cizalla. La figura 4.1 muestra la curva experimental y el correspondiente ajuste obtenidos por Steinwachs et al. [19] en este experimento, así como un pequeño diagrama que ilustra la configuración del reómetro utilizado. Asimismo, la figura 4.3 muestra las curvas resultantes del experimento para ambas concentraciones de fibras, recibidas de la Universidad KU Leuven. 23 4.1.2. Contracción vertical ante extensión horizontal El gel es depositado sobre la superficie activada con sulfo-SANPAH de un plato flexible de polidimetilsiloxano. Tras la polimerización, el gel se estira uniaxialmente en la dirección horizontal X a ritmo del 6 % por hora utilizando un motor paso a paso, y al mismo tiempo se mide con un microscopio la altura del gel (dirección Z), la cual es función del estiramiento aplicado. Durante todo el proceso se impide el movimiento en la dirección transversal Y. La expresión del tensor gradiente de deformación que caracteriza este experimento viene dada por, F=  λh0 0 0 1 0 0 0 λv (4.2) donde λhyλvrepresentan la extensión horizontal y la contracción vertical, respectivamente. La figura 4.2 muestra la curva experimental y el correspondiente ajuste obtenidos por Steinwachs et al. [19] en este experimento, así como un pequeño diagrama que ilustra la configuración utilizada. La figura 4.4 muestra la curva promedio de las resultantes para cada concentración de fibras (recibidas de la Universidad KU Leuven), puesto que éstas no presentan diferencias significativas entre sí. Figura 4.1: Curva experimental y correspondiente ajuste para el ensayo de cizalla simple (Steinwachs). Figura 4.2: Curva experimental y correspondiente ajuste para el ensayo de contracción vertical ante extensión horizontal (Steinwachs). 0 0.5 1 1.5 2 2.5 3 γ 0 20 40 60 80 100 120 σ12 [Pa] 0.6 mg/ml 1.2 mg/ml Figura 4.3: Curvas experimentales provenientes del ensayo de cizalla simple (KU Leuven). 0.95 1 1.05 1.1 1.15 1.2 1.25 1.3 λh 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 λv Figura 4.4: Curva experimental proveniente del experimento de contracción vertical ante extensión horizotal (KU Leuven). 24 4.2. Resultados En el caso del ensayo de cizalla simple, los ajustes se han llevado a cabo para el rango de deformaciones definido por γ∈[0,0.5] debido a que, al no presentar los modelos considerados puntos de inflexión en el rango de deformaciones del experimento, estos no son capaces de aproximar el comportamiento del gel más allá del punto en el que las curvas empíricas manifiestan un cambio de curvatura, alrededor de γ= 0.5. A partir de este punto, el comportamiento del gel entra en régimen plástico, manifestando una física diferente no contemplada por los modelos. Asimismo, rangos de deformación tan grandes no resultan de gran interés porque no son de aplicación en este problema. En el caso del ensayo de contracción vertical, el rango considerado es λh∈[1,1.05] debido a que es donde se obtienen mejores ajustes. 4.2.1. Neo-Hookeana Las figuras 4.5 y 4.7 muestran los mejores ajustes obtenidos en el caso del ensayo de cizalla simple para los geles de 0.6 y 1.2 mg/ml, respectivamente. Puede apreciarse claramente que el modelo no captura adecuadamente el comportamiento no lineal del gel. Esto se debe principalmente al carácter isocórico del gradiente de deformación en cizalla simple (J= 1), que oblitera toda contribución de κen el tensor de tensiones. De este modo, la tensión de cizalla adquiere la forma σ12 =µγ, comportándose entonces de forma lineal con respecto a la deformación de cizalla. El único parámetro ajustable a través de este ensayo es µ, que representa la pendiente de la recta resultante, y que adquiere sendos valores 20.1884 y26.7160 Pa. En las figuras 4.6 y 4.8 se representan los mejores ajustes para ambos geles en el caso del experimento de contracción vertical. Manteniendo los valores de µhallados a través del otro experimento se obtienen los valores de κ:1150 y1930 Pa, respectivamente. Se observa que el modelo sólo representa fielmente la porción de curva comprendida en el rango γ∈[1,1.025] aproximadamente. La curva teórica entonces se separa y, mientras que la curva experimental es cóncava en el rango aproximado γ∈[1,1.09], la anterior se presenta convexa desde el principio. En efecto, para valores de κmucho mayores que µ, la expresión (3.13) verifica que, si κµ−→ λv≈1 λh ,(4.3) siendo la curva teórica virtualmente independiente de los parámetros materiales y de curvatura estrictamente convexa. Es evidente que la curva proporcionada en este caso por el modelo Neo- Hookeano nunca podrá aproximar el comportamiento de la curva experimental más allá del punto de separación. 4.2.2. Holzapfel Las figuras 4.9 y 4.11 ilustran los ajustes obtenidos del ensayo de cizalla simple para ambos geles. Se aprecia un considerable aumento de precisión en la captura del comportamiento del gel con respecto al modelo Neo-Hookeano -incluso más notable en el de mayor concentración de fibras-, causado por la presencia de la función de tipo Fung (exponencial) en la definición del potencial elástico. Los valores obtenidos para las constantes características del modelo de Holzapfel son, respectivamente, k1= 16.5425 Pa y k2= 4.1789; y k1= 24.3827 Pa y k2= 3.5363. En cambio, ambos ajustes revelan un valor de µnulo, lo cual implica, teniendo en cuenta que κno tiene influencia en el comportamiento del gel ante cizalla simple, que sólo el conjunto de las fibras, y no la matriz, contribuye en la obtención de las tensiones. Esta observación puede tener fundamento en la hipótesis de deformación afín entre la red de fibras y la matriz planteada por Steinwachs et al. en la confección de su modelo [19], la cual es tanto más cierta cuanto mayor es la concentración de fibras del gel en cuestión. 25 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 γ 0 2 4 6 8 10 12 14 16 σ12 [Pa] Neo-Hook Experimental Figura 4.5: Ajuste para cizalla simple, 0.6 mg/ml - µ= 20.1884 Pa. 1 1.05 1.1 1.15 λh 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 λv Neo-Hook Experimental Figura 4.6: Ajuste para contracción vertical, 0.6 mg/ml - µ= 20.1884 Pa, κ= 1150 Pa. 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 γ 0 2 4 6 8 10 12 14 16 18 20 σ12 [Pa] Neo-Hook Experimental Figura 4.7: Ajuste para cizalla simple, 1.2 mg/ml - µ= 26.7160 Pa. 1 1.05 1.1 1.15 λh 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 λv Neo-Hook Experimental Figura 4.8: Ajuste para contracción vertical, 1.2 mg/ml - µ= 26.7160 Pa, κ= 1930 Pa. Si bien el modelo de Holzapfel se ajusta relativamente bien a las curvas del ensayo de cizalla simple, no ocurre lo mismo con el ensayo de contracción vertical, tal y como se desprende de las figuras 4.10 y 4.12. En concreto, manteniendo los valores de µ,k1yk2obtenidos previamente, y para κ= 970 Pa, el modelo produce curvas casi idénticas a las obtenidas con el modelo Neo-Hookeano. Esto sugiere que la contribución del término relativo a las fibras en la obtención de la contracción vertical es insignificante en comparación con aquella del relativo al parámetro κ. A este respecto, para valores de κsuficientemente grandes, la curva de contracción vertical es independiente de los parámetros materiales y opuesta en curvatura a la obtenida experimentalmente, de forma análoga a cómo ocurre en el caso del modelo Neo-Hookeano. Es por este motivo que las curvas son tan similares, y que el valor de κes el mismo para ambas concentraciones. En efecto, eliminando el término relativo a la matriz del potencial elástico y estudiando la curva de contracción vertical que produciría el mantener el término relativo a las fibras, se verifica la necesidad de conferir a k2valores extremadamente grandes a fin de producir una curvatura cóncava significativa en la gráfica, resultando k1irrelevante en la obtención de la misma (figura 4.13). El orden de magnitud de estos valores carece de sentido, pues trastocarían completamente los ajustes realizados para el ensayo de cizalla simple. 26 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 γ 0 2 4 6 8 10 12 14 16 18 σ12 [Pa] Holzapfel Experimental Figura 4.9: Ajuste para cizalla simple, 0.6 mg/ml - µ= 0 Pa, k1= 16.5425 Pa, k2= 4.1789. 1 1.05 1.1 1.15 λh 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 λv Holzapfel Experimental Figura 4.10: Ajuste para contracción vertical, 0.6 mg/ml - µ= 0 Pa, κ= 970 Pa, k1= 16.5425 Pa, k2= 4.1789. 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 γ 0 5 10 15 20 25 σ12 [Pa] Holzapfel Experimental Figura 4.11: Ajuste para cizalla simple, 1.2 mg/ml - µ= 0 Pa, k1= 24.3827 Pa, k2= 3.5363. 1 1.05 1.1 1.15 λh 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 λv Holzapfel Experimental Figura 4.12: Ajuste para contracción vertical, 1.2 mg/ml - µ= 0 Pa, κ= 970 Pa, k1= 24.3827 Pa, k2= 3.5363. 27 —˙ σ∗6=R˙ σRT, véase el capítulo cuarto de [20] para más información al respecto. Sin embargo, a la hora de linealizar (5.3), ABAQUS aproxima el término de fuerzas internas mediante la derivada de Jaumann o derivada co-rotacional, la cual es objetiva, DZV τ:δεdV [∆u] = ZV D[τ:δε][∆u]dV = =ZVh∇δu:[J c J− H ]:∇δu+τ:h[∇δu]T∇∆uiidV (5.4) en donde c Jes el tensor de elasticidad de Jaumann, que puede obtenerse como, J c J=J c +H, Hijkl =1 2[δjkτil +δikτjl +δjlτik +δilτjk](5.5) siendo c el tensor de elasticidad espacial o Cuarto Tensor de Elasticidad de Cauchy, c =1 Jφ∗[C], cijkl =1 JFiIFjJ FkK FlLCIJKL (5.6) en donde se ha empleado la notación de Einstein. A la hora de escribir la sub-rutina UMAT, hay que definir explícitamente el tensor de tensiones de Cauchy en función del tensor gradiente de deformación, y la matriz de rigidez tangente DDSDDE, la cual coincide con la cantidad J c J. Aprovechando las simetrías de ambas magnitudes, éstas habrán de ser expresadas en la versión particular de la notación de Voigt que ABAQUS emplea, σ=         σ11 σ22 σ33 σ12 σ13 σ23         DDSDDE =         D1111 D1122 D1133 D1112 D1113 D1123 D2211 D2222 D2233 D2212 D2213 D2223 D3311 D3322 D3333 D3312 D3313 D3323 D1211 D1222 D1233 D1212 D1213 D1223 D1311 D1322 D1333 D1312 D1313 D1323 D2311 D2322 D2333 D2312 D2313 D2323         (5.7) 5.2. Tensor de tensiones de Cauchy Recordando la formulación en configuración mixta de la ley constitutiva dada en la sección 2.2.3, el tensor de tensiones de Cauchy puede obtenerse aplicando la operación de empuje sobre el Primer Tensor de Piola-Kirchhoff, σ=1 Jφ∗[P] = 1 JP F T=1 J1 4πZ2π 0Zπ 0 Ψ0(λ)[F eΩ]⊗eΩ kF eΩksin θdθdφ·FT(5.8) aproximándose la integral sobre la superficie de la esfera unitaria mediante el método de cuadratura de Gauss-Legendre descrito en la sección 2.2.2. 5.3. Tensor de elasticidad material A diferencia con el resto, el tensor de elasticidad espacial no puede obtenerse directamente a través del potencial elástico. Alternativamente, dicho tensor ha de obtenerse a través de la relación que guarda con el tensor de elasticidad material (5.6). Recuperando la expresión (2.15), puede obtenerse el tensor de elasticidad en configuración mixta como, 34 A=∂P(λ(F),F) ∂F=∂P ∂λ ⊗∂λ ∂FΩ +∂P ∂F:41(5.9) que en notación indicial, AiJkL =∂PiJ ∂FkL =∂PiJ ∂λ ∂λ ∂FkL Ω +∂PiJ ∂Fmn δmkδnL = =eΩJeΩLkF eΩkΨ00(λ)−Ψ0(λ) kF eΩk2·[F eΩ]i[F eΩ]k kF eΩk+ Ψ0(λ)1 kF eΩkδikΩ (5.10) Puede probarse ([21]) que los tensores mixto y material se relacionan mediante, AiJkL =δikSJL +FiN CNJMLFkM (5.11) Haciendo uso de una de las simetrías menores de C, se tiene que, FiN CNJMLFkM =FiN CNJLM FkM −→ F·C·FT(5.12) No obstante, el tensor F·C·FTno posee, en general, simetrías menores, por lo que en el contexto de (5.11) es necesario reflejar dicha transposición, A=ξ+hF·C·FTik↔L,AiJkL =ξiJkL +hF·C·FTiiJLk , ξiJkL =δikSJL (5.13) donde la notación k↔Lindica transposición entre el par de índices kyL. Entonces, puede deducirse una expresión explícita para Ca través del despeje de la expresión anterior, C=F−1·[A−ξ]K↔l·[FT]−1,CIJKL =F−1 Ii [A−ξ]iJlK[FT]−1 lL (5.14) 5.3.1. Procedimiento alternativo Una forma alternativa, quizás más sencilla, de obtener el tensor de elasticidad material, es mediante la reformulación de la ley constitutiva en configuración material. Recordando la definición del cuarto invariante en la sección 2.2.2, puede extraerse que, I4=eT Ω·C·eΩ=eT Ω·FTF·eΩ=kF eΩk2(5.15) observación que permite reescribir λen términos del tensor de Cauchy por la derecha, λ(F) = kF eΩk−1−→ λ(C)=(I4)1/2−1(5.16) De este modo, la relación constitutiva en configuración material resulta, S= 2 ∂Ψ(λ(C)) ∂CΩ = 2 ∂Ψ ∂λ ∂λ ∂CΩ =[λ+ 1]−1·Ψ0(λ)·[eΩ⊗eΩ]Ω(5.17) Atendiendo a la expresión (2.9), uno puede obtener el tensor de elasticidad material como, C= 2∂S(λ(C)) ∂C= 2∂S ∂λ ⊗∂λ ∂CΩ (5.18) donde, ∂S ∂λ =−[λ+ 1]−2·Ψ0(λ)+[λ+ 1]−1·Ψ00(λ)·[eΩ⊗eΩ]Ω(5.19) ∂λ ∂C=1 2[λ+ 1]−1·[eΩ⊗eΩ](5.20) 35 5.4. Validación de la implementación En esta sección se presentan algunas pruebas realizadas con el fin de comprobar la validez de la implementación. Las pruebas en concreto han sido extraídas de la guía proporcionada por ABAQUS para la programación de sub-rutinas UMAT [22]. Pueden distinguirse dos tipos esenciales de pruebas: aquellas diseñadas para verificar la implementación de la ley constitutiva, y aquellas diseñadas para verificar la implementación del Jacobiano, o matriz de rigidez tangente. Esta pruebas se efectúan sobre un cubo unitario mallado con un único elemento. Los resultados obtenidos tras la realización de las pruebas son entonces comparados con los obtenidos a través de la implementación de la ley constitutiva en MATLAB, a fin de verificar que la sub-rutina funciona correctamente. Los parámetros materiales empleados en dichas pruebas son arbitrarios y elegidos por conveniencia, en concreto, el parámetro de pandeo dose ha elegido con un valor lo suficientemente grande como para evitar el plateau característico de la curva de compresión que produce el modelo de Steinwachs, y así solventar potenciales problemas de convergencia. El cojunto de valores empleado es: (κo, do, λs, ds) = (200 Pa,1,0.01,5). 5.4.1. Verificación del algoritmo de integración de las tensiones En estas pruebas se definen todos los desplazamientos existentes en la deformación en la sección de condiciones de contorno. De este modo, la simulación no depende de la definición del Jacobiano ya que el tensor gradiente de deformación está determinado desde el principio, y la obtención de las tensiones es directa. Extensión uniaxial Se parte del tensor gradiente de deformación, F=  1.500 0 1 0 0 0 1  (5.21) que proporciona el siguiente tensor de tensiones de Cauchy, σ=  24.1106 0 0 0 4.1838 0 0 0 4.1838  [Pa](5.22) Los resultados de la simulación revelan un estado tensional homogéneo idéntico a (5.22), tal y como ilustra la figura 5.1 en el caso de la componente σ11. Aunque por compacidad no se presenten el resto de figuras, se ha comprobado la coincidencia de las seis componentes de la tensión de Cauchy. Extensión uniaxial en dirección oblicua Se parte del tensor gradiente de deformación, F=  1.5 0 0 0.125 1 0 0.25 0 1  (5.23) que caracteriza un estado de extensión uniaxial en la dirección [√21/8]·[0.5,0.125,0.25]T. El tensor de tensiones de Cauchy asociado es, σ=  24.8369 3.0366 6.0732 3.0366 4.6281 0.6150 6.0732 0.6150 5.5507  [Pa](5.24) 36 Los resultados de la simulación revelan un estado tensional homogéneo idéntico a (5.24), tal y como ilustra la figura 5.2 en el caso de la componente σ13. Aunque por compacidad no se presenten el resto de figuras, se ha comprobado la coincidencia de las seis componentes de la tensión de Cauchy. (Avg: 75%) S, S11 +2.411e+01 +2.411e+01 +2.411e+01 +2.411e+01 +2.411e+01 +2.411e+01 +2.411e+01 +2.411e+01 +2.411e+01 +2.411e+01 +2.411e+01 +2.411e+01 +2.411e+01 Step: Step−1 Increment 10: Step Time = 1.000 Primary Var: S, S11 Deformed Var: U Deformation Scale Factor: +1.000e+00 ODB: UMAT.odb Abaqus/Standard 6.13−3 Fri Sep 04 21:36:30 Hora de verano romance 2020 X Y Z Figura 5.1: σ11 en extensión uniaxial. (Avg: 75%) S, S13 +6.073e+00 +6.073e+00 +6.073e+00 +6.073e+00 +6.073e+00 +6.073e+00 +6.073e+00 +6.073e+00 +6.073e+00 +6.073e+00 +6.073e+00 +6.073e+00 +6.073e+00 Step: Step−1 Increment 10: Step Time = 1.000 Primary Var: S, S13 Deformed Var: U Deformation Scale Factor: +1.000e+00 ODB: UMAT.odb Abaqus/Standard 6.13−3 Fri Sep 04 22:25:28 Hora de verano romance 2020 X Y Z Figura 5.2: σ13 en extensión uniaxial en dirección oblicua. 5.4.2. Verificación de la precisión del Jacobiano En este caso, se define un estado de tensiones en la sección de condiciones de contorno con el fin de obtener los desplazamientos correspondientes. Tracción-Compresión Se parte del siguiente estado tensional, 37 σ=  7.1545 0 0 0−3.6356 0 0 0 −2.8580  [Pa](5.25) asociado al tensor gradiente de deformación, F=  1.2 0 0 0 0.85 0 0 0 0.9 (5.26) que provoca una tracción en la dirección X, y dos compresiones en las direcciones YyZ. Los desplazamientos de tres puntos cualesquiera en cada una de las tres caras frontales recibirán los valores, uX=  0.2 0 0 ;uY=  0 −0.15 0 ;uZ=  0 0 −0.1 (5.27) Los resultados de la simulación revelan un campo de desplazamientos homogéneo idéntico a (5.27), tal y como ilustra la figura 5.3 en el caso de la componente u2(dirección Y). Aunque por compacidad no se presenten el resto de figuras, se ha comprobado la coincidencia de los desplazamientos en las otras dos direcciones. U, U2 −1.501e−01 −1.376e−01 −1.251e−01 −1.126e−01 −1.000e−01 −8.754e−02 −7.504e−02 −6.253e−02 −5.002e−02 −3.752e−02 −2.501e−02 −1.251e−02 +0.000e+00 Step: Step−1 Increment 10: Step Time = 1.000 Primary Var: U, U2 Deformed Var: U Deformation Scale Factor: +1.000e+00 ODB: UMAT.odb Abaqus/Standard 6.13−3 Fri Sep 04 22:10:36 Hora de verano romance 2020 X Y Z Figura 5.3: Desplazamiento en la dirección Yen un estado de tracción-compresión. 38 Capítulo 6 Conclusión Este Trabajo Fin de Máster se ha desarrollado en el seno de una colaboración con investigadores de la Universidad KU Leuven, una oportunidad única para contribuir con el trabajo realizado, aunque humildemente, a un proyecto de gran envergadura: la caracterización mecánica de hidrogeles de colágeno que aproximen con precisión las especificidades de la llamada matriz extracelular, entorno hospedador de la célula y factor determinante en su evolución; un pilar fundamental en el estudio de la técnica TFM. Desde un punto de vista personal, la consecución de este trabajo ha supuesto un notable complemento a los conocimientos adquiridos durante el transcurso de mis estudios universitarios, debido a que ha requerido el estudio profundo de conceptos avanzados de la mecánica del sólido deformable, desde meras intuiciones físicas hasta desarrollos matemáticos complejos y abstractos. En esta sección final se expone una recapitulación de los principales resultados y observaciones presentados en el trabajo, así como las futuras tareas que siguen al mismo. 6.1. Resumen de los resultados Se han ajustado con éxito los datos experimentales provenientes de los ensayos de cizalla simple y contracción vertical frente a extensión horizontal —efecto Poisson—, obtenidos por colaboradores en la Universidad de Lovaina, empleando tres leyes de comportamiento hiperelásticas: Neo-Hook, Holzapfel y Steinwachs. De los resultados se extrae que, aparentemente, la mejor opción para la caracterización de los hidrogeles es el modelo de Steinwachs. El modelo Neo-Hookeano falla a la hora de capturar el comportamiento no-lineal de los hidrogeles al mostrar una dependencia lineal entre la tensión y la deformación de cizalla. Si bien el modelo de Holzapfel solventa este problema —ajusta con precisión dicha curva—, sólo lo hace parcialmente, ya que falla por completo, al igual que la ley Neo-Hookeana, al predecir la curva de contracción vertical. Se ha verificado que la influencia del término de Holzapfel en la obtención de la última es insignificante en comparación con la del término Neo-Hookeano, por lo que ambos modelos producen curvas casi idénticas. El modelo de Steinwachs es el único que logra capturar con relativa precisión la curva de contracción vertical, a la vez que ajusta fielmente la curva de cizalla simple. No obstante, este último presenta algunos problemas que dificultan su implementación en problemas de mayor envergadura. El modelo de Steinwachs se cimenta sobre la asunción de que en un hidrogel la deformación de fibras y matriz se produce de forma independiente. Bajo esta hipótesis, el potencial elástico del hidrogel se concibe como el promedio de los potenciales elásticos de todas las fibras, excluyéndose la contribución de la matriz acuosa. En este sentido, debido a que la rigidez de una fibra a compresión se establece prácticamente nula —la fibra pandea—, el comportamiento promedio que ofrece el modelo muestra también rigidez nula ante compresión. La rigidez nula a compresión es fuente de abundantes problemas de carácter computacional: convergencia, distorsión de malla, pérdida de convexidad, entre otros. Asimismo, la hipótesis de deformaciones afines sólo es válida si la concentración de 39 fibras es lo suficientemente grande. En caso de que la concentración de fibras sea pequeña, la rigidez general del hidrogel será comparable a la de la matriz, y entonces sería necesario considerar la contribución de la misma en el comportamiento general del material. 6.2. Trabajos futuros Los siguientes pasos a desarrollar en el marco de esta colaboración son: Corroboración de la participación de la matriz en la deformación de los hidrogeles, sirviéndose de datos correspondientes a ensayos de compresión que se espera recibir de Lovaina en el futuro próximo. Además, se espera recibir información acerca de la deformación de las fibras en el interior del gel, para así poder verificar la hipótesis de deformaciones afines. Desarrollo de un modelo que caracterice mejor el comportamiento del término de la fibra, a partir de los modelos de Holzapfel y Steinwachs, y que también incluya la contribución de la matriz en la definición de la función densidad de energía de deformación. En consecuencia, se espera poder conseguir mejores ajustes de las curvas experimentales en un rango más amplio de deformaciones, incluyendo la región de compresión. Este nuevo modelo sería por sí mismo una contribución interesante en el campo de la TFM. Una vez que el modelo se encuentre definido, y comprobadas las mejoras que conlleva, se adaptará la sub-rutina UMAT ya escrita para que incluya dicha ley constitutiva. Utilización de la subrutina desarrollada para su uso en problemas reales de TFM, en colaboración con la Universidad KU Leuven, para la investigación del microentorno mecánico de células ante diversas condiciones fisiológicas. 40 Bibliografía [1] F. Spill, D. S. Reynolds, R. D. Kamm, M. H. Zaman. Impact of the physical microenvironment on tumor progression and metastasis. Current Opinion in Biotechnology, Volumen 40, Páginas 41–18. Elsevier. Agosto de 2018. DOI: https://doi.org/10.1016/j.copbio.2016.02.007 [2] F. Spill, C. Bakal, M. Mak. Mechanical and Systems Biology of Cancer. Computational and Structural Biotechnology Journal, Volumen 16, Páginas 237–245. Elsevier. 2018. DOI: https: //doi.org/10.1016/j.csbj.2018.07.002 [3] U. S. Schwarz, J. R. D. Soiné. Traction force microscopy on soft elastic substrates: A guide to recent computational advances. Biochimica et Biophysica Acta (BBA) - Molecular Cell Research, Volumen 1853, Edición 11, Parte B, Páginas 3095–3104. Elsevier. Noviembre de 2015. DOI: https://doi.org/10.1016/j.bbamcr.2015.05.028 [4] S. A. Maskarinec, C. Franck, D. A. Tirrell, G. Ravichandran. Quantifying cellular traction forces in three dimensions. Proceedings of the National Academy of Sciences of the United States of America, Volumen 106, Número 52, Páginas 22108–22113. National Academy of Sciences. Diciembre de 2009. DOI: https://doi.org/10.1073/pnas.0904565106 [5] L. G. Griffith, M. A. Swartz. Capturing complex 3D tissue physiology in vitro. Nature Reviews - Molecular Cell Biology, Volumen 7, Páginas 211–224. Nature. Marzo de 2006. DOI: https: //doi.org/10.1038/nrm1858 [6] J. Toyjanova, E. Bar-Kochba, C. López-Fagundo, J. Reichner, D. Hoffman-Kim, C. Franck. High Resolution, Large Deformation 3D Traction Force Microscopy. PLoS ONE 9(4). Abril de 2014. DOI: https://doi.org/10.1371/journal.pone.0090976 [7] Y. J. Suh. Cancer Cell Mechanics. Biofluidics Lab. Cornell University. Web. Enlace: https: //biofluidics.bee.cornell.edu/cancer-cell-mechanics.html [8] Y. Huang, C. Schell, T. B. Huber, A. N. Şimşek, N. Hersch, R. Merkel, G. Gompper, B. Sabass. Traction force microscopy with optimized regularization and automated Bayesian parameter selection for comparing cells. Scientific Reports, Volumen 9, Número 539. Nature. Enero de 2019. DOI: https://doi.org/10.1038/s41598-018-36896-x [9] P. C. Hansen. Rank-Deficient and Discrete Ill-Posed Problems: Numerical Aspects of Linear Inversion. Mathematical Modeling and Computation. Society for Industrial and Applied Mathematics. 1998. DOI: https://doi.org/10.1137/1.9780898719697 [10] G. A. Holzapfel. Nonlinear Solid Mechanics: A Continuum Approach for Engineering. John Wiley & Sons. 2000. [11] J. Bonet, R. D. Wood. Nonlinear Continuum Mechanics For Finite Element Analysis. Cambridge University Press. 1997. [12] R. W. Ogden. Non-Linear Elastic Deformations. Dover Publications. Mineola, Nueva York. 1997. 41 [13] T. J. Pence, K. Gou. On compressible versions of the incompressible neo-Hookean material. Mathematics and Mechanics of Solids, Volumen 20, Número 2, Páginas 157–182. SAGE journals. Febrero de 2015. DOI: https://doi.org/10.1177/1081286514544258 [14] B. A. Lane, K. A. Harmon, R. L. Goodwin, M. J. Yost, T. Shazly, J. F. Eberth. Constitutive modeling of compressible type-I collagen hydrogels. Medical Engineering & Physics, Volumen 53, Páginas 39–48. Elsevier. Marzo de 2018. DOI: https://doi.org/10.1016/j.medengphy. 2018.01.003 [15] O. Gültekin, H. Dal, G. A. Holzapfel. On the quasi-incompressible finite element analysis of anisotropic hyperelastic materials. Computational Mechanics, Volumen 63, Páginas 443–453. Nature. Marzo de 2019. DOI: https://doi.org/10.1007/s00466-018-1602-9 [16] G. A. Holzapfel, H. W. Weizsäcker. Biomechanical behavior of the arterial wall and its numerical characterization. Computers in Biology and Medicine, Volumen 28, Edición 4, Páginas 377–392. Elsevier. Julio de 1998. DOI: https://doi.org/10.1016/S0010-4825(98)00022-5 [17] Y. C. Fung, K. Fronek, P. Patitucci. Pseudoelasticity of arteries and the choice of its mathematical expression. American Journal of Physiology — Heart and Circulatory Physiology, Volumen 237, Número 5, Páginas H620–H631. American Physiological Society. Noviembre de 1979. DOI: https://doi.org/10.1152/ajpheart.1979.237.5.H620 [18] K. Atkinson. Numerical Integration On The Sphere. Journal of the Australian Mathematical Society, Volumen 23, Edición 3, Páginas 332–347. Cambridge University Press. Enero de 1982. DOI: https://doi.org/10.1017/S0334270000000278 [19] J. Steinwachs, C. Metzner, K. Skodzek, N. Lang, I. Thievessen, C. Mark, S. Münster, K. E. Aifantis, B. Fabry. Three-dimensional force microscopy of cells in biopolymer networks. Nature Methods, Volumen 13, Páginas 171–176. Nature. Diciembre de 2015. DOI: https: //doi.org/10.1038/nmeth.3685 [20] E. W. V. Chaves. Notes on Continuum Mechanics. Engineering Mechanics Lecture Notes on Numerical Methods in Engineering and Sciences. Springer. 2013. [21] "Nonlinear finite elements/Rate form of hyperelastic laws". Wikiversity. Wikimedia Foundation. Web. 2017. Enlace: https://en.wikiversity.org/wiki/Nonlinear_finite_ elements/Rate_form_of_hyperelastic_laws [22] ABAQUS Online Guide. Lecture 6: Writing a UMAT or VUMAT. Writing User Subroutines with ABAQUS. Simulia ABAQUS. Web. Enlace: https://imechanica.org/files/ Writing%20a%20UMAT.pdf [23] O. C. Zienkiewicz, R. L. Taylor, J. Z. Zhu. The Finite Element Method: Its Basis and Fundamentals, séptima edición. Butterworth-Heinemann, Elsevier. 2013. [24] Greg von Winckel (2020). Legendre-Gauss Quadrature Weights and Nodes, MATLAB Central File Exchange. Accedido en septiembre de 2020. Enlace: https://www.mathworks.com/ matlabcentral/fileexchange/4540-legendre-gauss-quadrature-weights-and-nodes [25] Mathworks (2020). lsqnonlin: Solve nonlinear least-squares (nonlinear data-fitting) problems (r2016a). Accedido en septiembre de 2020. Enlace: https://www.mathworks.com/help/ optim/ug/lsqnonlin.html 42