Acerca de la precisión en la determinación de las constantes de difusión en hidrogeles mediante simulación numérica

Armando Blanco ^1,2*, Gema González ¹, Euro Casanova ², María Pirela ^1,3*, Alexander Briceño ³

¹ Centro de Ingeniería de Materiales y Nanotecnología. Instituto Venezolano de Investigaciones Científicas. e-mail: gemagonz@ivic.gob.ve, mpirela@ivic.gob.ve

² Departamento de Mecánica. Universidad Simón Bolívar. e-mail: ajblanco@usb.ve, ecasanov@usb.ve

³ Centro de Química. Instituto Venezolano de Investigaciones Científicas. e-mail: abriceno@ivic.gob.ve

RESUMEN

La precisión de los esquemas numéricos empleados para discretizar la ecuación de difusión, es de gran importancia en la asignación de la constante de difusividad que mejor representa la penetración de un solvente en un hidrogel. Se comparan los valores obtenidos para la difusividad de un solvente en un hidrogel, determinada mediante dos metodologías diferentes, desarrolladas utilizando esquemas basados en diferencias finitas y elementos finitos, basadas en la discretización de la ecuación de difusión que modela el proceso de penetración del solvente en el hidrogel, determinando, al mismo tiempo, el cambio de volumen del mismo. La aplicación se realiza en muestras de hidrogeles de Poliacrilamida. Las constantes de difusividad fueron obtenidas por ajuste de las predicciones numéricas con los resultados experimentales. Los resultados alcanzados permiten concluir que importantes variaciones en la constante de difusión, en una misma muestra, se obtuvieron dependiendo de la metodología numérica utilizada. En consecuencia, una revisión profunda de las metodologías utilizadas para estimar los distintos parámetros asociados con los modelos de difusión debe ser realizada. Adicionalmente, se analizó la poca influencia de las incertidumbres de las distintas variables en la estimación de la difusividad, los valores específicos de la relación radio/espesor que corresponden al proceso de hinchamiento más lento y el efecto de las variaciones en la concentración de equilibrio en la dinámica del hinchamiento.

Palabras clave: Hinchamiento de hidrogeles, Diferencias finitas, Elementos finitos, Coeficiente de difusividad.

On the accuracy in the determination of diffusivity constants in hydrogels by numerical simulation

ABSTRACT

The accuracy of the numerical schemes used for discretizing the diffusion equation is of great importance for determining the value of the diffusivity constant that best represents the diffusion of a solvent in a hydrogel. In this paper we compare the values obtained for the diffusivity constant of a solvent in a hydrogel, by using two different methodologies based on finite difference and finite element schemes respectively. These methodologies consider the discretization of the diffusion equation to model the process of penetration of the solvent in the hydrogel, while modeling the volume change of the hydrogel. The specific application of both approaches is performed on cylindrical samples of polyacrylamide hydrogels. Diffusivity constants were determined by fitting the numerical predictions with experimental results. Numerical results show that huge differences in the values of diffusivity constant were obtained depending on the numerical method used, requiring, therefore, a thorough review of the methodologies used so far. Additionally, little influence of the uncertainty of parameters in diffusivity law was observed in the estimated diffusivity; specific values of the radius / thickness ratio corresponding to slower swelling process were calculated, and finally the effect of variations in the equilibrium concentration in swelling dynamics was determined.

Keywords: Swelling of hydrogels, Finite difference, Finite elements, Diffusivity coefficients.

Recibido: agosto 2013 Recibido en forma final revisado: marzo 2014

INTRODUCCIÓN

Los hidrogeles son materiales poliméricos hidrofílicos, conformados por una red de cadenas entrecruzadas constituidas por macromoléculas entrelazadas con segmentos de homopolímeros y copolímeros hidrofílicos (Sotoudeh et al. 2010) capaces de absorber cantidades importantes de agua sin disolverse (Kou et al. 1990). Un hidrogel completamente seco, esto es, sin ninguna traza de agua en su interior, se denomina Xerogel (Pal et al. 2009). Cuando los hidrogeles están hinchados, son blandos, elásticos, presentan alta afinidad por el agua, alta estabilidad térmica y mecánica así como biocompatibilidad, aspectos que los hace atractivos para aplicaciones en diversas áreas tales como la ingeniería química, medicina, ingeniería biomédica, farmacia, alimentos y agricultura (Sotoudeh et al. 2010).

El proceso de hinchamiento de un hidrogel ha sido descrito por varios autores (Ganji et al. 2010). Una vez que la muestra entra en contacto con el agua, ésta comienza a penetrar en su interior. Durante este proceso, el agua ocupa espacios entre las cadenas poliméricas, separándolas. Esta separación permite a su vez que nuevas cantidades de agua penetren y ocupen más espacios. Sin embargo, debido al entrecruzamiento de las cadenas, y a la elasticidad de las mismas, se alcanza un estado de equilibrio, en el cual la presión osmótica ejercida por el agua se balancea con la resistencia elástica de la estructura (Ganji et al. 2010). Durante este proceso de absorción de agua o hinchamiento, los hidrogeles incrementan su volumen hasta alcanzar su máxima capacidad de almacenamiento o saturación de agua (Kou et al. 1990).

La dinámica del proceso de hinchamiento depende de múltiples factores tales como las características particulares de los materiales constituyentes del hidrogel y la forma específica de la muestra.

Para el modelado del hinchamiento de hidrogeles, básicamente dos tipos distintos de familias de modelos matemáticos han sido ampliamente utilizados. Una primera familia está constituida por modelos matemáticos, que sólo consideran el comportamiento global de la muestra de hidrogel en términos de la cantidad de solvente, generalmente agua, absorbida. Este tipo de modelos considera tanto la difusión del líquido como la relajación de la estructura polimérica de acuerdo con la expresión (Peppas & Sahlin, 1989):

donde: M_t y M_∞ son las masas de agua en el instante de tiempo t y en condiciones de equilibrio (t infinito) respectivamente, k₁ y k₂ son constantes asociadas con los procesos de difusión y relajación de la estructura del polímero respectivamente, cuyos valores dependen del sistema particular solventepolímero y n es denominado como el “exponente de descarga” (en inglés release exponent) y caracteriza la dinámica del proceso de hinchamiento. Debido a la forma particular de la expresión utilizada para exponer la dependencia funcional de la Ecuación (1), este modelo se denomina “modelo de ley de potencias”. Dependiendo de los valores que tome n, el proceso es catalogado como difusión “Fickeana” (n =1/2), Caso II (n =1), absorción anómala (1/2 < n <1) (De Kee et al. 2005). Rara vez se caracterizan procesos de hinchamiento con n < ½ (Ganji et al. 2010).

Entre las limitaciones de este tipo de modelos, se encuentra la imposibilidad de considerar los efectos debido a la forma particular de la muestra de hidrogel. Una segunda familia de modelos que permite sobrepasar estas limitaciones, se basa en las leyes de difusión de Fick. En estos modelos, la concentración de solvente en la matriz polimérica es modelada siguiendo la ecuación de difusión clásica (Crank, 1979):

donde: C es la fracción volumétrica local de solvente y D es la difusividad. Resultados excelentes han sido obtenidos cuando la difusividad es representada utilizando una dependencia exponencial en la concentración de solvente, como la propuesta por Fujita (1961). Dado que la difusividad no es constante y debido a la expansión de las muestras de hidrogel durante el proceso de hinchamiento, no es posible encontrar soluciones analíticas a la ecuación de difusión, en estas condiciones, lo que obliga a considerar el uso de técnicas numéricas para su resolución. En particular, excelentes resultados con aplicaciones numéricas basadas en el método de las diferencias finitas, que consideran muestras cuyas formas corresponden a discos cilíndricos, han sido obtenidos por autores como Siepmann et al. (1999).

Sin embargo, aplicaciones basadas en esquemas de discretización, cuya base es el método de las diferencias finitas, tienen como limitación principal la forma de las muestras que pueden ser consideradas, restringiéndose, en general, a formas tipo disco para las cuales un dominio de cálculo simplificado puede definirse.

Recientemente, un esquema de discretización para la ecuación de difusión, basado en el método de los elementos finitos, que considera la determinación precisa de la ubicación de los nodos de la malla para cada instante del proceso de hinchamiento, ha sido planteado por Blanco et al. (2013), para aplicaciones específicas en el hinchamiento de hidrogeles. Este esquema fue probado con excelentes resultados para muestras de formas distintas tales como discos cilíndricos, cápsulas y conjuntos cilindro-cono.

Sin embargo, hasta donde es del conocimiento de los autores, la comparación de las predicciones de la dinámica del hinchamiento de hidrogeles, entre modelos basados en esquemas en diferencias finitas y elementos finitos, no ha sido realizada. En este artículo se realiza una comparación entre ambas familias de métodos. La forma de las muestras consideradas corresponde a discos cilíndricos, debido a que son las únicas que pueden ser representadas en mallas para las cuales es posible utilizar fácilmente esquemas en diferencias finitas. La comparación con experimentos realizados en laboratorio, permitieron calibrar y comparar ambos esquemas, en especial, en lo que corresponde a las estimaciones de los coeficientes de difusividad. Adicionalmente, la influencia de la morfología de las muestras (relación diámetro/espesor) en el proceso de hinchamiento, así como el efecto de la incertidumbre en la determinación de las constantes adimensionales en la expresión global del coeficiente de difusión, es analizada utilizando el modelo basado en elementos finitos.

En este artículo en la sección 2 se formula el modelo matemático utilizado para modelar el hinchamiento de hidrogeles. Luego, en la siguiente sección se discuten brevemente la aplicación de ambos esquemas, diferencias finitas y elementos finitos, a la ecuación de difusión con dominio móvil y coeficientes de difusividad variables. Además, se presentan tres esquemas distintos de generación de mallas, los cuales fueron utilizados con ambos modelos para representar el dominio de cálculo en las distintas fases del proceso de hinchamiento, siendo su principal aporte la validación de ambos modelos así como el conjunto de los casos modelados.

MODELO MATEMÁTICO

El modelado matemático de la penetración de un solvente en una matriz polimérica inicialmente seca, se realizó considerando que el proceso puede ser representado por la Ecuación de difusión (2), que escrita en coordenadas cilíndricas se expresa como:

donde: C = C(r,θ,z,t) es la concentración de solvente en el punto (r,θ, z) en el instante de tiempo t; r,θ y z son las coordenadas radial, azimutal y axial respectivamente, y D=D(C) es la difusividad del solvente en la matriz polimérica. La difusividad D, al igual que en Siepmann et al. (1999) y Blanco et al. (2013), es considerada como en Fujita (1961), expresada como:

donde: D_eq y C_eq son el coeficiente de difusión y la concentración del solvente una vez que el estado de hinchamiento de equilibrio es alcanzado respectivamente, y β es una constante adimensional que caracteriza la dependencia de la difusividad con la concentración de D.

Consideremos que la muestra de hidrogel completamente seca, cuya forma corresponde a un disco cilíndrico de radio r₀ y espesor h₀, es definida tal como muestra la Figura 1.

Dada la simetría de la geometría de la muestra, la concentración C(r,θ,z,t) = C(r,z,t) es independiente de la variable θ y, en consecuencia, la Ecuación (3) puede ser escrita como:

La Ecuación (5) debe ser resuelta con la difusividad expresada por la Ecuación (4), con las condiciones iniciales y de borde apropiados. Inicialmente, en t=0, puesto que la muestra está completamente seca, la concentración inicial de solvente en el interior corresponde a C(r,z,0)=0. Así mismo, en las fronteras r= r₀ y z= ±h_0/2, C(r₀,z,0) = C_eq y C(r,h₀,0) = C_eq, respectivamente.

La imposición de las condiciones de frontera requiere que si se define la frontera Ω como la superficie del disco cilíndrico, entonces C(Ω(t),t)=C_eq. Esta condición será expresada de manera distinta dependiendo de la estrategia utilizada para deformar el dominio físico a medida de que se produce el hinchamiento.

MODELOS NUMÉRICOS

Como ya ha sido mencionado, en este trabajo se comparan dos modelos distintos basados en los métodos de diferencias finitas y elementos finitos. A continuación se describen brevemente ambos modelos.

Modelo en diferencias finitas

Basados en la propuesta original de Siepmann et al. (1999), la Ecuación (3) es discretizada en diferencias finitas de primer orden en el tiempo y segundo orden en el espacio. El esquema escogido es de tipo explícito. A partir del conocimiento de los valores de la concentración del solvente en el instante de tiempo t, los nuevos valores en el instante t+Δt pueden ser calculados directamente. Una vez conocidas las concentraciones en el instante t+Δt, puede determinarse la cantidad de masa total de solvente que ingresó a la matriz polimérica y, en consecuencia, el nuevo volumen del conjunto gel-solvente puede ser calculado. Se considera, al igual que Siepmann et al. (1999), que el volumen de agua que ingresa a la matriz polimérica por cada una de las direcciones radial y axial conduce a incrementos de volumen en esas mismas direcciones, de manera de que el incremento de volumen en cada dirección es proporcional al área superficial respectiva. Luego, se calculan los nuevos valores del radio y el espesor del disco y se genera una nueva malla, conservando el mismo número de nodos en cada dirección.

Dadas las condiciones de simetría propias del disco cilíndrico considerado, el dominio de cálculo puede ser reducido a un rectángulo, tal como se muestra en la Figura 2.

Modelo en elementos finitos

El modelo numérico basado en el método de los elementos finitos fue presentado por Blanco et al. (2013) y se describe brevemente para efectos de la integridad del artículo.

La formulación en elementos finitos se realizó utilizando el método de los residuos ponderados (Zienkiewicz et al. 2005). El residuo es definido a partir de la Ecuación (5) como:

Este residuo es minimizado imponiendo:

donde: ψ es una función de ponderación arbitraria, la cual debe anularse en fronteras en las que se imponen condiciones de Dirichlet y el elemento de volumen, dv=2πrdA, dado que el dominio de cálculo es axisimétrico.

Introduciendo la Ecuación (6) en la (7) y después de simplificar las integrales la formulación débil del problema es obtenida:

Utilizando el enfoque de Galerkin (Zienkiewicz et al. 2005) la concentración y las funciones de peso son interpoladas a partir de sus valores en los nodos utilizando las mismas funciones matriciales de interpolación:

Con estas interpolaciones, la Ecuación (8) se transforma para quedar como:

donde: M^e, K^e y f_(t) ^e son matrices y vectores elementales:

donde: B representa la matriz de las derivadas de la función de interpolación.

Finalmente, la formulación en elementos finitos del problema es obtenida ensamblando todos los elementos en una única ecuación:

Mċ_(t) + Kc_(t) = f_(t)                  (14)

En este trabajo, el modelo en elementos finitos fue construido utilizando elementos triangulares lineales axisimétricos de tres nodos y un grado de libertad por nodo, correspondiente a la concentración en ese punto. Este tipo de elementos fue escogido, por una parte, debido a la simplicidad de la implementación en el modelo numérico y, por otra parte, dado que las propias incertidumbres asociadas con la determinación de los valores de las distintas variables que describen el proceso de hinchamiento (contenido de agua, forma de las muestras) no permiten diferenciar entre las soluciones numéricas procedentes de elementos finitos distintos.

Generación del dominio de cálculo de las mallas

Para deformar el dominio de cálculo, se consideraron dos metodologías distintas. En primer lugar, se utilizó la misma metodología que para el modelo en diferencias finitas, descrita en el apartado anterior. Sin embargo, el cálculo se realizó empleando elementos finitos triangulares, tal como se muestra en la Figura 3.

En segundo lugar, se utilizó la misma metodología propuesta por Blanco et al. (2013). Ésta consiste en generar un sistema de ecuaciones algebraicas no lineales el cual es resuelto numéricamente en cada iteración. Cada nodo de la malla de cálculo es referenciado por los índices i y j. Las consideraciones que generarán el sistema de ecuaciones se obtienen de las cuatro condiciones siguientes:

(a) Los nodos se desplazan a lo largo de líneas directrices fijas que corresponden a valores de pendiente constante en el tiempo, tal como se muestran en la Figura 4.

(b) El volumen de cada elemento se calcula en cada instante de tiempo en función de la cantidad de solvente que ingresa a ese elemento particular.

(c) Las líneas correspondientes a valores j=cte son perpendiculares al eje de simetría z.

(d) Dado que en este trabajo sólo se compararon ambos modelos para discos cilíndricos, las líneas correspondientes a valores j=cte son perpendiculares al eje de simetría radial.

La Figura 4 presenta mallas iníciales típicas (t=0) y deformadas (t>0). La Figura 4b muestra como a medida de que el solvente es absorbido por el hidrogel, los elementos de la malla se deforman dependiendo de la concentración local del solvente.

Verificación de los modelos numéricos

Dado que no existe solución analítica para validar ambos modelos numéricos, con la finalidad de verificar la precisión de los mismos, se realizaron dos simulaciones, en dominios

computacionales fijos, de la difusión de solvente en un disco típico, considerando que el coeficiente de difusión era (a) constante y (b) definido de acuerdo con la Ecuación (4). El ingreso de solvente en un disco cilíndrico de radio 4,88mm y espesor 8,93mm fue calculado suponiendo que el coeficiente de difusión se expresaba como:(a) D_eq=5,2.10^-2 mm²/min, β=0 y C_eq=0,9135 y (b) D_eq=5,2.10^-2mm²/min, β=2,5 y C_eq=0,9135. En particular, el valor de β es el mismo utilizado por Siepman et al. (1999) y la sensibilidad de los distintos modelos numéricos ante variaciones del mismo, es considerada en la sección de resultados. Estos valores corresponden a una muestra de las utilizadas por Blanco et al. (2013). Se analizó tanto las predicciones numéricas para la absorción global de solvente (agua) por el disco, como los tiempos de llegada del agua a cuatro puntos distintos del disco, en cada caso.

La Figura 5, presenta la evolución de la masa relativa de agua absorbida por el disco en función del tiempo, construida a partir de la ecuación reportada para medir el índice de hinchamiento o porcentaje de hidratación S (Katime et al., 2004), 100*(Masa(t=t)-Masa(t=0))/ Masa(t=0), para ambos casos. Resulta claro (Figura 5) que las diferencias entre las predicciones numéricas no son significativas.

La Figura 6 presenta la localización de cuatro puntos en el cuerpo del disco para comparar la concentración de agua en función del tiempo para ambos modelos.

Las Figuras 7 y 8, presentan la evolución local de la concentración del solvente en los cuatro puntos seleccionados, considerando el coeficiente de difusión constante (Figura 7) y variable (Figura 8).

Tanto en la Figura 7 como en la 8 se aprecian tendencias análogas. En particular, se observa que los tiempos predichos por ambos modelos numéricos para la llegada de agua son bastante similares en cada uno de los puntos considerados.

En consecuencia, dado que tanto la descripción global del hinchamiento como la determinación de la concentración local de agua en distintos puntos del disco cilíndrico son bastante cercanas, se considera que ambos esquemas numéricos han sido verificados.

PREPARACIÓN DE LAS MUESTRAS

Para comparar las predicciones de los esquemas en diferencias finitas y elementos finitos, se utilizó un conjunto de discos cilíndricos preparados con los siguientes materiales y métodos. Los hidrogeles de Poliacrilamida (PAAm), obtenidos por polimerización en solución vía radical libre de monómeros de Acrilamida (AAm, 100 % m/m) utilizando N, N´, Metilenbisacrilamida (MBAAm, 1,1 % m/m) como agente entrecruzante y en presencia de Persulfato de Amonio (PSA, 2 % m/m) como iniciador de la reacción (iniciación química). El tiempo de disolución, en 0,7 ml de agua destilada, para la AAm y MBAAm fue de 5 minutos cada uno y luego se adicionó el PSA. Las síntesis se realizaron en un baño termostatado, utilizando un agitador magnético con calefacción, a una temperatura de 30±1°C y bajo agitación continua de 800 rpm, hasta la formación de los hidrogeles (monolitos). Estos monolitos posteriormente, fueron hidratados durante cuatro (4) días con la intención de lavarlos y así extraer los monómeros no reaccionantes durante las síntesis; luego cortados en discos de diferentes tamaños y puestos en una estufa a 40°C, para así llevarlos al estado xerogel. Las curvas de hinchamiento fueron obtenidas de acuerdo con la metodología descrita por Katime & Rodríguez, (2001) y Rojas de Gascue et al. (2010).

COMPARACIÓN DE LOS MODELOS NUMÉRICOS

Determinación de los coeficientes de difusividad

Los modelos basados en diferencias finitas y elementos finitos fueron ajustados considerando las curvas de hinchamiento obtenidas de un conjunto de discos preparados según los procedimientos descritos en la sección anterior. La Tabla 1, presenta las características de cada muestra.

El ajuste de los modelos numéricos fue realizado únicamente sobre el valor del coeficiente de difusión. El valor de β se fijó en 2,5 mientras que la concentración en equilibrio, C_eq, se determinó a partir del hinchamiento máximo alcanzado por cada muestra (Katime et al. 2004).

La Figura 9 presenta las curvas de hinchamiento para ocho de las muestras consideradas. Es notorio que ambos modelos, una vez que la constante de difusión ha sido ajustada, son capaces de reproducir de manera bastante precisa el proceso global de captura de solvente.

Sin embargo, desde un punto de vista cualitativo, el modelo de elementos finitos basado en generación de malla radial describe mejor el hinchamiento. Esto es confirmado al comparar los coeficientes de correlación entre la data experimental y las predicciones numéricas presentadas en la Tabla 2.

Así mismo, la Tabla 2 presenta los coeficientes de difusión obtenido para cada muestra. Es de hacer notar la gran diferencia en los valores obtenidos para este coeficiente (a) entre las predicciones numéricas de ambos esquemas (Tabla 2, columna 2 con columnas 3 y 4) y (b) entre las obtenidas con el modelo de elementos finitos y ambos tipos de mallas (Tabla 2, columnas 3 y 4 entre sí). Dado que ambos esquemas fueron extensamente verificados considerando una malla fija, y que las mismas también aparecen utilizando el mismo esquema numérico (Tabla 2, columnas 3 y 4 entre sí), estas diferencias sólo pueden ser ocasionadas por las especificidades del proceso de generación de los dominios computacionales y la generación de las nuevas mallas en cada instante de tiempo.

Puesto que la generación de mallas considerando los cambios locales de masa de solvente es menos restrictiva, al no imponer condiciones particulares a la forma del hidrogel durante el proceso de hinchamiento, permitiendo que la muestra se deforme libremente, es de esperar que las predicciones numéricas utilizando esta metodología sean más precisas que las obtenidas a partir de mallas que condicionen la deformación de las muestras.

Este postulado encuentra soporte adicional cuando se observa que, en efecto, durante el proceso de hinchamiento, la forma del hidrogel cambia apreciablemente como se observa en la Figura 10, para una de las muestras utilizadas en este estudio.

Claramente, durante las etapas iniciales del disco no se hincha proporcionalmente, al contrario, el mayor crecimiento se realiza en la zona de intersección entre las caras planas y la periferia circular, dado que es allí donde se produce un mayor influjo de agua. Sin embargo, en el estado de total hinchamiento el disco recupera su forma original.

Resulta de especial importancia destacar las desviaciones encontradas entre las predicciones numéricas de los coeficientes de difusión entre distintos esquemas, aún cuando los procesos de hinchamiento puedan reproducirse con alta precisión. Esta desviación podría ser de gran importancia en procesos en los cuales se requiera de alta precisión en la captura del solvente por parte del hidrogel. Las causas de esta desviación podrían deberse a imprecisiones en la dinámica de generación de las mallas, o a fallas en el modelo matemático en la descripción del proceso de deformación durante los instantes iniciales de hinchamiento al no considerar los efectos del propio movimiento de la malla.

Análisis de influencia de la relación de aspecto del disco

Dependiendo de la forma del disco cilíndrico, esto es, del valor de la relación radio/espesor, distintos comportamientos en lo que respecta a la dinámica de la absorción de solvente pueden esperarse. Así, en discos muy achatados (radio/ espesor>>1) o muy esbeltos (radio/espesor<<1), el solvente puede alcanzar rápidamente todos los espacios del interior.

Para relaciones de forma intermedias el proceso de captura de solvente debe ser en principio más lento. En particular, debe existir un valor de la relación radio/espesor para la cual este proceso de absorción se realiza de la manera más lenta posible.

Para determinar este valor de la relación de aspecto y su dependencia con el tamaño del disco, se modeló el proceso de toma de solvente considerando dos discos cilíndricos base, de masa fija y valores típicos de las propiedades de las muestras antes analizadas. La Tabla 3, presenta las propiedades de cada disco.

La relación de aspecto fue modificada para cada disco, conservando la cantidad de masa constante. Luego, definiendo un valor para el volumen del disco, el valor del radio del disco seco era determinado al fijar el del espesor (o viceversa). La Figura 11, presenta las curvas de hinchamiento en función del tiempo para cada caso base con relaciones de aspecto que oscilan entre 0,1 y 2,7.

De la Figura 11 se desprende, tal como se espera, el hinchamiento se produce más rápidamente en muestras achatadas o esbeltas. Sin embargo, aún cuando el proceso de hinchamiento se realiza más rápidamente en la muestra de menor masa, en ambos casos, la relación de aspecto para la absorción más lenta estuvo entre 0,3375 y 0,8, independientemente de la cantidad de masa de la muestra. En consecuencia, para muestras cuya relación de aspecto se encuentre entre los valores arriba mencionados, es de esperar que la utilización de un modelo como el descrito por la Ecuación (1) permita ajustar de manera bastante precisa el proceso de hinchamiento global en el tiempo.

Análisis de influencia de la precisión en la determinación de β

En general, las constantes de difusividad y la constante adimensional β se determinan a través del ajuste estadístico de curvas de hinchamiento obtenidas por medios experimentales. Con la finalidad de determinar el efecto de imprecisiones en el procedimiento de ajuste para la determinación del valor de la constante adimensional β, se realizó un conjunto de simulaciones del hinchamiento de una disco cilíndrico (muestra 1, Tabla 3) considerando variaciones de β de ±10% y ±20% en torno al valor utilizado comúnmente para agua, β =2,5 (Siepmann et al. 1999). La dinámica del hinchamiento en función del tiempo, para los cinco casos considerados, se presenta en la Figura 12.

De la Figura 12 se desprende que incluso desviaciones en los valores de β de orden 20% tienen poco efecto global en el proceso de hinchamiento del hidrogel. En particular, en el caso considerado la máxima desviación respecto al valor promedio del porcentaje de hinchamiento, %S, que es de 6,67%.

Análisis de influencia de Ceq en la dinámica del hinchamiento

La dinámica del proceso de hinchamiento está condicionada, de manera importante, por el valor de la concentración de equilibrio del solvente, C_eq. Mientras mayor es esta constante, mayor es la capacidad de captación de agua por parte del hidrogel. Sin embargo, debido a la no linealidad introducida por la difusividad, la dinámica del hinchamiento del hidrogel no es lineal. Los efectos del valor de C_eq se estudiaron para la muestra 1 de la Tabla 3, variando los valores de C entre 0,71 y 0,84. La Figura 13, presenta las curvas de hinchamiento en función del tiempo y normalizadas en función del máximo valor de hinchamiento posible para cada caso.

En la Figura 13a se observa, como era de esperarse, una mayor velocidad de hinchamiento a medida de que el valor de C_eq es mayor. Sin embargo, en términos relativos al valor máximo del hinchamiento, S_max, mostrados en la Figura 13b, el proceso se desarrolla más rápidamente para muestras cuyos valores de C_eq es menor. En consecuencia, si se desea recoger una cantidad dada de solvente rápidamente, podría ser mejor estrategia introducir en lugar de un disco con un valor alto de C_eq, una mayor cantidad de discos con valores más bajos de C_eq.

CONCLUSIONES

En este trabajo se compararon las estimaciones en los coeficientes de difusión que caracterizan el hinchamiento de muestras de hidrogeles, cuya forma correspondía a discos cilíndricos, utilizando dos modelos numéricos distintos basados en los métodos de las diferencias finitas y de los elementos finitos, así como dos esquemas de generación de mallas.

Ambos modelos numéricos son capaces de reproducir con gran precisión el proceso de captura del solvente por parte del hidrogel. No obstante, tanto desde un punto de vista cualitativo como cuantitativo, el ajuste del modelo basado en elementos finitos con deformación arbitraria de la muestra fue más preciso que el de diferencias finitas.

Diferencias importantes en las predicciones numéricas fueron encontradas dependiendo de la metodología utilizada para reproducir la forma del hidrogel durante el proceso del hinchamiento, aún cuando se emplea el mismo modelo numérico, en este caso específico, el método de los elementos finitos. En particular, al considerar los valores obtenidos para el coeficiente de difusión utilizando los dos esquemas de discretización y ambas metodologías de generación de mallas, se observaron diferencias en las predicciones del coeficiente de difusión que alcanzaban, para algunas muestras, casi el 100% respecto al valor promedio obtenido, tal como se observa en la Tabla 2. Adicionalmente, las predicciones de los coeficientes de difusión entre ambos métodos numéricos difieren hasta en un 46% respecto al valor medio obtenido utilizando ambos modelos numéricos. Estas diferencias sólo pueden ser debidas a las metodologías utilizadas para representar la deformación del hidrogel y para generar la malla de cálculo.

Además se demostró que imprecisiones en el ajuste de la constante β, en la expresión de la difusividad tienen poca importancia práctica en la predicción de la dinámica del proceso de hinchamiento. Así mismo, se comprobó que variaciones en los valores de la concentración de equilibrio no conducen a cambios lineales en la dinámica del hinchamiento e inclusive se mostró que la dinámica de hinchamiento es relativamente más rápida en aquellos casos en los cuales la concentración de equilibrio es más baja.

En futuros trabajos se analizarán en detalle tanto los efectos de los métodos utilizados para representar la deformación del hidrogel y desplazar los nodos de la malla, como la incertidumbre en la determinación experimental de las constantes de difusividad debido a parámetros tales como la ubicación de la muestra en el monolito y pequeñas variaciones en las condiciones de realización de los experimentos.

REFERENCIAS

1. Blanco, A., González, G., Casanova, E., Pirela, M., Briceño, A. (2013). Mathematical Modeling of Hydrogels Swelling Based on the Finite Element Method, Applied Mathematics, 4(8), 2013; pp. 161-170. doi: 10.4236/am.2013.48A022.        [ Links ]

2. Crank, J. (1979). The Mathematic of Diffusion. 2da Ed. Clarendon Press, Oxford.        [ Links ]

3. De Kee, D., Liu, Q., Hinestroza, J. (2005). Viscoelastic (Non-Fickian) Diffusion. The Canadian Journal of Chemical Engineering, 83(6); pp. 913-929.        [ Links ]

4. Fujita, H. (1961). Diffusion in polymer-diluent systems. In Fortschritte der Hochpolymeren-Forschung, Springer Berlin Heidelberg; pp. 1-47.        [ Links ]

5. Ganji, F., Vasheghani-Farahani, S., Vasheghani-Farahani, E. (2010). Theoretical description of hydrogel swelling: a review. Iranian Polymer Journal, 19(5); pp. 375-398.        [ Links ]

6. Katime, I., Katime, T. O., Katime, T. D. (2004). Los materiales inteligentes de este milenio: Hidrogeles macromoleculares. Síntesis, propiedades y aplicaciones. Editorial de la Universidad del País Vasco. Bilbao. España.        [ Links ]

7. Katime, I. & Rodríguez, E. (2001). J. Macromol. Sci. Pure Appl. Chem., A38 (5 & 6): 543-558.        [ Links ]

8. Kou, J. H., Fleisher, D., Amidon, G. L. (1990). Modeling drug release from dynamically swelling poly (hydroxyethyl methacrylate-co-methacrylic acid) hydrogels. Journal of Controlled Release, 12(3); pp. 241-250.        [ Links ]

9. Pal, K., Banthia, A. K., Majumdar, D. K. (2009). Polymeric hydrogels: characterization and biomedical applications. Designed monomers and polymers, 12(3), 197-220.        [ Links ]

10. Peppas, N. A. & Sahlin, J. J. (1989). A simple equation for the description of solute release. III. Coupling of diffusion and relaxation. International Journal of Pharmaceutics, 57(2); pp. 169-172.        [ Links ]

11. Rojas de Gáscue, B., Ramírez, M., Prin, J.L., Torres, C., Bejarano, L., Villarroel, H., Rojas, L., Murillo, M., Katime, I. (2010). Hidrogeles de Acrilamida/ Ácido Acrílico y de Acrilamida/Poli(Ácido Acrílico): Estudio de su capacidad de Remediación de Efluentes industriales. Revista Latinoamericana de Metalurgia y Materiales, 30 (1): 28-39.        [ Links ]

12. Siepmann, J., Podual, K., Sriwongjanya, M., Peppas, N. A., Bodmeier, R. (1999). A new model describing the swelling and drug release kinetics from hydroxypropyl methylcellulose tablets. Journal of pharmaceutical sciences, 88(1); pp. 65-72.        [ Links ]

13. Sotoudeh, S., Pourfallah, G., Barati, A., Davarnejad, R., Farahani, M. A., Memar, A. (2010). Dynamical Modeling and Experimental Analysis on the Swelling Behavior of the sIPN Hydrogels. Industrial & Engineering Chemistry Research, 49(20); pp. 10111-10115.        [ Links ]

14. Zienkiewicz, O. C., Taylor, R. L., Zhu, J.Z. (2005). The Finite Element Method: Its Basis and Fundamentals. Capítulo 3, 6ta Ed. Oxford (Inglaterra). Elsevier.        [ Links ]