El propósito principal de la presente publicación consiste en la determinación adecuada del orden fraccional en el modelo viscoelástico de Zener y el análisis de las implicaciones que se derivan de la precisión en la obtención del mismo. Esto se realiza empleando el método numérico de Levenberg-Marquardt a partir de valores reportados en la literatura. Los parámetros utilizados son obtenidos para un segmento de arteria empleando un experimento de relajación a los esfuerzos. Una vez determinado el orden fraccional se procede a determinar la solución del modelo empleando la función de dos parámetros de Mittag-Leffler y la operación de convolución, con la finalidad de comparar el comportamiento del modelo de orden fraccional con el de orden entero e identificar sus principales diferencias. Se obtienen además las deformaciones que experimenta la arteria bajo el estímulo de un pulso sanguíneo normal y otros disímiles, simulando el efecto provocado en el flujo sanguíneo por una arritmia y el ocasionado por el proceso de ventilación mecánica. Por último se analiza la respuesta dinámica del material bajo una serie de pulsos, utilizando la operación de convolución y el método numérico de Gauss-Kronrod, identificando la precisión en la reproducción de los mismos en el modelo de orden fraccional comparado con el de orden entero.
The main purpose of this publication is the proper determination of the fractional order at the Zener viscoelastic model's and the analysis of the implications derived from the accuracy in obtaining this. The procedure is performed using the numerical method of Levenberg-Marquardt from values reported in the literature. The parameters used are obtained from an artery segment using a stress relaxation test. After determining the fractional order is proceeded to find the model solution using, the function of Mittag-Leffler with two parameters and the convolution operation in order to compare the behavior of the fractional model vs the integer order and identify their key differences. The displacements present on the artery are obtained, under the stimulus of a normal blood pulse and two dissimilar, simulating the effect in blood flow caused by an arrhythmia and by mechanical ventilation process. Finally the dynamic response from the pulses is analyzed using the convolution operation and Gauss Kronrod numeric method, where the solution's accuracy obtained by the fractional model is observed, an compared with the integer order model.
Determinar los fenómenos mecánicos que ocurren en las arterias resulta de suma importancia debido a su influencia en la fisiología arterial y en el tratamiento-seguimiento de los principales padecimientos arteriales ocasionados por los efectos del flujo sanguíneo [1]. Estos padecimientos se deben generalmente a: la calcificación de las paredes arteriales, los aneurismas [2], la obstrucción de las arterias por acumulación de placas de colesterol, padecimiento comúnmente conocido como aterosclerosis, [3] y por los efectos perjudiciales de la hipertensión [4], por citar algunos de los principales.
Los materiales que forman el cuerpo humano se dividen principalmente en dos grandes ramas, los tejidos blandos y los mineralizados. En la primera rama se encuentran las arterias, venas, músculos, piel, tendones, entre otros. La segunda se compone principalmente de los huesos y los dientes. Aunque en algunos casos suelen presentarse, tejidos blandos con formaciones calcificadas, estos no son habituales.
Las arterias a su vez se componen principalmente de redes de colágeno y elastina, además de tejido muscular liso lo que les proporciona un tipo de comportamiento material sumamente complejo.
Generalmente en la mayoría de los estudios biomecánicos de los tejidos blandos, estos han sido caracterizados bajo los modelos de comportamiento material viscoelástico [5]. Debido a que poseen características propias de los materiales elásticos, generalmente sólidos, y a su vez también se comportan como un fluido viscoso.
Una de las principales características de los materiales elásticos consiste en su capacidad de almacenar energía cuando son deformados bajo cargas, y en devolver toda esta energía cuando la carga se deja de aplicar. De forma contraria en el fluido viscoso la energía mecánica es continuamente disipada en forma de calor. Sin embargo existen materiales, como el caso de este estudio, que simultáneamente almacenan y disipan energía mecánica cuando se sujetan a un proceso de cargas, este tipo de materiales se conocen como viscoelásticos. Por lo que en estos no solamente se relaciona la tensión mecánica con la deformación, sino que además involucra su razón de cambio con respecto al tiempo.
Los materiales viscoelásticos suelen ser modelados a través de arreglos consistentes de resortes y amortiguadores, como lo son los modelos reológicos de Maxwell y de Kelvin-Voigt. Los modelos mencionados anteriormente son propuestos inicialmente por Fung [6] para su aplicación en la modelación de tejidos blandos. Sin embargo sólo se consideraban como aproximaciones teóricas, ya que estos modelos tradicionales presentan únicamente una caracterización parcial o descriptiva del fenómeno. Los mismos han sido modificados usualmente a partir de los modelos clásicos agregando más elementos en serie o paralelo al original, logrando una mejor caracterización del comportamiento mecánico.
El problema derivado de realizar esta adición de elementos es que se aumenta considerablemente la complejidad del modelo haciendo poco conveniente su uso.
Sin embargo en los últimos años un antiguo concepto con novedosas aplicaciones, el cálculo fraccional, ha sido utilizado en la modelación y caracterización de materiales viscoelásticos, principalmente los polímeros.
Al parecer la viscoelasticidad es el campo de mayor aplicación de los operadores diferenciales e integrales fraccionales... se demuestra que el uso de derivadas fraccionales para el modelado matemático de materiales viscoelásticos resulta natural... Podlubny [7]
También existen una serie de investigaciones orientadas a la caracterización de los efectos mecánicos en los tejidos biológicos blandos empleando modelos fraccionales viscoelásticos como lo son: el utilizado en tejido extraído del hígado [8], la determinación del comportamiento mecánico del tejido graso empleando núcleos derivativos [9] o el del tejido cerebral [10] a partir de los modelos clásicos agregando un nuevo elemento fraccional conocido como spring-pot, cuya derivada es de orden no entero y se encuentra para dichos modelos entre cero y uno [11].
A su vez los modelos viscoelásticos de orden fraccional son empleados en la obtención de los cambios de las propiedades mecánicas de tejidos blandos para la detección temprana de diversas situaciones patológicas, principalmente el cáncer de próstata [12], donde el papel de la determinación del orden fraccional en el modelo juega un rol sumamente importante debido a que la exactitud del mismo está directamente asociado con la certeza del procedimiento. Por lo que en diversas investigaciones se desarrollan metodologías en función de su adecuada obtención [13,14].
2Cálculo fraccionalEl origen del cálculo fraccional data desde el mismo fundamento del cálculo diferencial [7] cuando en una carta enviada a Leibniz el 30 de Septiembre de 1695 en la que L’Höpital le cuestiona el significado que tendría su notación desarrollada para la n-ésima derivada dnydxn, el qué sucedería en el caso n=12, a lo que Leibniz le responde:
“Usted puede ver por eso señor que uno puede expresar por una serie infinita una cantidad como d1/2ydx1/2. Aunque las series infinitas y geométricas son relaciones distantes, las series infinitas admiten sólo el uso de exponentes que son enteros, no hace todavía el uso de exponentes fraccionarios... esto conduciría a una paradoja, de la que algún día se extraerán consecuencias útiles”
2.1La derivada fraccional definida en ℝ+El operador diferencial local de la derivada de orden n para una variable independiente t, Dtn=dndtn es solamente el operador inverso por la izquierda del operador integral no local del n-campo integral aItn teniendo como punto de inicio cualquier valor a<t. De hecho, para cualquier función continua con derivada continua, es decir, f(t)∈ C2 se reconoce:
ydonde f(k)(a+) es la k-ésima derivada de la función f evaluada en a por la derecha. Como consecuencia de esto, tomando a=0 se requiere que 0Dtα sea definida como inversa por la izquierda de 0Itα. Para dicho propósito primero se introduce el número entero:y entonces se define la derivada fraccional de Riemann-Liouville de orden α>0:como:para complementar la definición de la derivada se tiene que 0Dt0=I. Donde Γ(x) es la función Gamma la cual se emplea como una generalización del factorial para el caso donde x no es un número entero.
Intercambiando en la ecuación (3) el proceso de derivación e integración, se llega a la conocida definición de la derivada fraccional de Caputo, tratada de forma extensa en [15], de orden α>0 definida como:
donde:
Para α∉ℤ la definición de Caputo, ecuación (5), se requiere que la m-ésima derivada de f(t) sea absolutamente integrable y no solamente por partes.
2.2La transformada de Laplace para derivadas fraccionalesCabe señalar que resulta de mayor utilidad la definición de derivada fraccional de Caputo en aquellos problemas propios de la física y la ingeniería donde se utilizan condiciones iniciales las cuales se encuentran expresadas en términos de derivadas de orden entero. Esto puede ser visualizado fácilmente mediante el uso de la transformada de Laplace, misma que requiere en su definición de transformada de una derivada, de las condiciones iniciales las cuales están dadas en términos de derivadas de orden entero. De hecho, análogamente a la definición para la transformada de una derivada de orden entero, para la derivada de Caputo de orden α con m−1<α≤m, se tiene:
dondepara la cual L;s es el operador de Laplace, y f˜(s) es la transformación de la función f(t) del dominio del tiempo al de la frecuencia. Sin embargo cuando las condiciones iniciales son cero, la ecuación (7), se reduce a:
La transformada de Laplace puede ser expresada en términos de las funciones de Mittag-Leffler, en particular
dondecon α,β∈ℝ+ y λ∈ℝ. Eα,β(−λtα) es la función de dos parámetros de Mittag-Leffler, la cual surge como una generalización de la función exponencial.3Modelo viscoelástico de Zener
El modelo viscoelástico de Zener consiste en una adecuación fraccional del sólido lineal estándar (SLS por sus siglas en inglés) [16]. Para lograr obtener la versión fraccional es necesario que se reemplace la derivada de primer orden por la derivada fraccional de orden α∈(0, 1) en su ecuación constitutiva. De aquí se obtienen las siguientes relaciones de esfuerzo deformación:
donde σ(t) es el esfuerzo, ϵ(t) es la deformación y a1, b1, m son constantes del modelo relacionadas con los coeficientes de elasticidad y viscosidad.
La teoría de la viscoelasticidad establece una relación entre el esfuerzo y la deformación a través del módulo de relajación a los esfuerzos G(t) y el de fluencia J(t) de la siguiente manera:
ydonde ξ es una variable de convolución. En el caso del modelo de Zener el módulo de fluencia tiene la forma:y el módulo de relajación a los esfuerzos se define como:además sustituyendo en la ecuación (12), es posible obtener el módulo complejo en el dominio de la frecuenciadonde E*(s) es el módulo complejo y se compone de las contribuciones de la parte elástica y la parte viscosa a la deformación.4Resultados
Utilizando la definición de la transformada de Laplace para derivadas fraccionales, ecuación (7) y la de la función de dos parámetros de Mittag-Leffler, ecuación (9), se obtiene:
donde λ=mb1.
Desarrollando los tres primeros términos de las sumatorias para la definición de las funciones de Mittag-Leffler, ecuación (10), se obtiene:
ydonde δ(t) es la función delta de Dirac y la introducción de la misma permite evadir la singularidad ocasionada por el segundo parámetro en la función de Mittag-Leffler. Aplicando la transformada de Laplace para la operación de convolución,
Es posible obtener una solución analítica de la deformación experimentada, bajo la aplicación de una función de tensión temporal [17].
Sea ϵs(t) la respuesta temporal correspondiente a una función escalón unitario σs(t) utilizada como la fuerza aplicada de manera temporal a la que está sometida la arteria, es posible obtener para el modelo de Zener, ecuación (16), la respuesta de deformación. La misma se muestra en la figura 1 para los valores de los parámetros a1=5.4871, b1=5.8712, m=0.68 y α=0.408.
Como es de esperarse, la solución del modelo de Zener exhibe el clásico comportamiento material que presentan los tejidos blandos, al ser estos materiales viscoelásticos. Donde la deformación en relación con la aplicación de la carga se presenta de manera no lineal mostrando una tendencia hacia los materiales considerados como fluidos viscosos. Sin embargo en la segunda etapa cuando la fuerza cesa en su aplicación, se muestra un estado de relajación, el cual es propio de los materiales viscoelásticos. Comparado con el modelo SLS el cual presenta a partir del momento donde se aplica inicialmente la carga un comportamiento elástico clásico es decir casi lineal y un tiempo de relajación más breve que el presentado por el modelo de Zener. Estas características del modelo de Zener, resultan sumamente favorables dado que el mismo se utilizará para simular los efectos mecánicos en un tejido blando, particularmente un segmento de arteria.
Los parámetros mecánicos que se utilizan en la presente investigación son los mismos que los obtenidos por Craiem et al [18] empleando un experimento de relajación a los esfuerzos, este se realiza sobre segmentos de arteria extraídos de cuatro pacientes. El mismo consiste en aplicar una carga cíclica periódica sobre el segmento de arteria y al cesar su aplicación se miden los tiempos que tarda en regresar a su estado tensorial original.
Se elige al primer paciente en el nivel de esfuerzo mecánico bajo debido a que la geometría de la arteria con la que se cuenta se construyó a partir de imágenes médicas obtenidas mediante un TAC para un paciente de 37 años de edad y estos parámetros son los que más se aproximan a las condiciones reales.
Como señalan los autores, estos no encontraron variaciones significativas entre los valores de las constantes mecánicas para los pacientes entre los dos modelos que emplearon, sin embargo si identificaron variaciones entre los valores de α[19], este último valor representa el orden de la derivada fraccional. Por esta razón se decidió realizar un ajuste sobre el valor de α, para determinar si este resultado es aquel que mejor se aproxima al conjunto de datos o si se debía considerar un nuevo valor del orden fraccionario y que implicaciones se derivan de realizar dicho ajuste.
La determinación del orden fraccional se realiza mediante la implementación de un algoritmo basado en el método numérico de Levenberg-Marquardt para ajuste de curvas mediante mínimos cuadrados para modelos no lineales, aproximando la función de relajación a los esfuerzos G(t) a un conjunto de puntos experimentales obtenidos en [19].
Los resultados obtenidos se muestran en el tabla (1) donde es posible observar que incrementando el número de datos utilizados en el método numérico se obtienen valores cercanos a 0.4, mismo que se encuentra dentro de los rangos adecuados para tejidos blandos como se muestra en diversos estudios, por ejemplo los de Craiem et al [20] que oscilan alrededor de 0.2 para el caso de un segmento de arteria, pero más cercanos a los presentados por Nagehan y Ergin [21] para otros tejidos blandos. Este valor de α indica que el material tiene características predominantes del tipo elástico al estar por debajo de 0.5, pero con una etapa de fluido viscoso al identificarse claramente el fenómeno de relajación a los esfuerzos, cuando la carga es retirada.
Relación entre el número de puntos de control utilizados al momento de implementar el método numérico, el valor que se obtiene para α y la calidad de la aproximación en términos de la desviación estándar de la curva obtenida contra la experimental.
| Cantidad de datos utilizados | Valor obtenido de α | Desviación estándar |
|---|---|---|
| 450 | 0.2132 | 0.0048 |
| 900 | 0.2336 | 0.0034 |
| 1,800 | 0.2542 | 0.0015 |
| 4,500 | 0.2813 | 0.0004 |
| 4,500,000 | 0.4026 | 0.0001 |
| 45,000,000 | 0.4081 | 0.00003 |
A continuación en la figura 2 se muestra un gráfico en el cual se presentan los distintos valores de α para el ajuste del modelo de Zener, resulta conveniente mencionar que a partir del valor α=0.4081 la aproximación aumenta su desviación estándar, es decir comienza a alejarse del mínimo absoluto en la regresión. Además se observa que a pesar de que el rango en el que se mueve la función es relativamente pequeño (0.72 ∼ 0.77), el ajuste propuesto inicialmente en la investigación de la cual se obtuvieron los parámetros dista considerablemente del nuevo valor de α propuesto. Y al estar trabajando con un sistema dinámico, es necesario tener en consideración el como se afecta la precisión del mismo, ya que el sistema puede tener un respuesta sensible a las condiciones iniciales, como se muestra adelante en este estudio.
Por otra parte el módulo de fluencia J(t) es comúnmente utilizado en experimentos donde se mantiene de forma constante el esfuerzo y se registran los diferenciales de la deformación en función del tiempo.
Los materiales viscoelásticos presentan un comportamiento variado, el mismo debe ser cuadrático en pequeños intervalos de tiempo y lineal en tiempos mayores, esta función se muestra en la figura 3. Aquí se puede apreciar inicialmente el comportamiento lineal que proporciona el modelo de orden entero y la notable diferencia entre los valores de α donde se puede apreciar como para el valor de α=0.2542 el módulo de fluencia entra rápidamente en el estado lineal lo cual no sucede para el valor de α=0.4081.
Como se mencionó anteriormente, el módulo complejo es el vector obtenido como suma de las contribuciones de los módulos elástico (parte real) y viscoso (parte imaginaria), en la figura 4 se muestra el módulo complejo en el dominio de la frecuencia s. Una vez más se puede apreciar la diferencia entre la respuesta obtenida para el modelo de orden entero y los de orden fraccional donde para estos últimos también es posible apreciar una considerable diferencia para bajas frecuencias.
Es importante mencionar que es a partir del módulo complejo que se obtienen los parámetros en el dominio de la frecuencia, los cuales son introducidos en un software especializado para el cálculo mediante elementos finitos, para la obtención de la distribución de esfuerzos en el segmento de arteria.
Por último se desea determinar cuál es la respuesta temporal de deformación para la arteria a partir de la aplicación de una serie de pulsos sanguíneos tal como sucede en la realidad. Primeramente se construye una función que asemeje el pulso con las fases de sístole (contracción) y diástole (relajación), descripciones del flujo arterial pueden ser consultadas en [22].
La función que simula el efecto del pulso sanguineo se describe a continuación:
y la forma del pulso se muestra en la figura 5.
Para este pulso unitario se obtiene la respuesta temporal del segmento de arteria, el mismo se muestra en la figura 6. Aquí es posible apreciar la precisión que presenta el modelo fraccional en comparación con el del SLS, en el tiempo inicial se observa el efecto del pulso con sus etapas para el modelo de Zener y en cambio para el SLS esta etapa se refleja en dos estados casi lineales, para después caer en un estado de relajación nuevamente lineal. A diferencia del fraccional donde existe una caída no lineal y un estado pronunciado de relajación.
Sin embargo es de sumo interés determinar los efectos del pulso como una sucesión de los mismos por lo que se genera una función periódica pπ(t), la cual se muestra en la figura 7 para los tres primeros períodos.
Aplicando el método numérico de Gauss-Kronrod en la solución de las integrales presentes en el proceso de convolución es posible obtener la respuesta temporal de la aplicación del tren de impulsos, la misma que se muestra en la figura 8 para su estado inicial, y para tiempos prolongados en la figura 9. Se puede observar en principio una zona de inestabilidad, pero alrededor de los cinco segundos se presenta una zona de oscilación que muestra básicamente la relación de tensión y deformación de interés en la investigación. Se debe tener en cuenta, como se mencionó anteriormente, que debido al uso de la definición de la transformada de Laplace para la derivada fraccional de Caputo, se requiere que las condiciones iniciales del sistema sean definidas para órdenes enteros de las derivadas. Con la finalidad de simplificar el sistema, las condiciones iniciales se toman como nulas. De aquí se concluye que la dinámica del sistema entra en un estado estable y este permanecería así de no retirar el estímulo sobre el mismo, figura 10.
Ahora el segmento de arteria será sometido a una carga “inusual” simulando un proceso de arritmia [23,24] el cual se compone de un pulso inicial en el que antes de haber terminado su período se presenta un nuevo pulso lo cual genera un traslape entre los dos. Este se define como paπ(t) donde los pulsos tienen una misma frecuencia pero distinta amplitud. La simulación del fenómeno se muestra en la figura 11
A diferencia de la estabilidad mostrada en el caso anterior, para el presente tren de impulsos con diferente amplitud se tiene una zona de inestabilidad que tiende aumentar a lo largo del tiempo, esto se muestra en la figura 12.
Inicialmente se muestra una zona de inestabilidad, la cual sobrepasa el estado de deformaciones del pulso regular en menos de tres segundos, sin embargo a partir del tiempo t=5 el sistema debería comenzar a estabilizarse, cosa que no sucede. Después del tiempo inicial continuá un estado creciente donde al cabo de 20 segundos prácticamente se duplica la magnitud de la deformación. Es claro que este fenómeno fisiológicamente no es posible, debido a que de presentarse esta condición por un período de tiempo similar se hubiese producido un infarto. Ver figura 13.
Realizando un acercamiento a la zona de inestabilidad señalada en la figura 13, se genera la figura 14 donde se puede observar el comportamiento repetitivo de la respuesta temporal al flujo que simula el proceso de arritmia y nuevamente es posible visualizar la sensibilidad que presenta el modelo fraccional al reproducir este efecto.
Para concluir se analiza el cambio en el flujo común provocado por un dispositivo de ventilación mecánica descrito en [25] donde se detectan picos de presión en conjunto con un desplazamiento oscilatorio de la señal en la aorta. La representación del flujo sanguíneo de la presión se describe en la figura 15.
La respuesta del sistema es similar a la presentada por el efecto de arritmia donde al igual que esta última se tiene una respuesta temporal que sobrepasa los límites presentados por el flujo normal de 120/80 mmHg y donde también se observa un efecto no transitorio de oscilaciones con una pendiente creciente. Para fines de la investigación este efecto se prolonga en un tiempo t=40s, sin embargo en la realidad este fenómeno sólo se presenta durante breves intervalos de tiempo en el proceso de ventilación mecánica. La respuesta a este comportamiento se describe en la figura 16.
Realizando un acercamiento a la zona de interés, figura 17, se observa nuevamente el comportamiento creciente y oscilatorio que se presenta en la respuesta del sistema dinámico, donde también se logra apreciar la sensibilidad que presenta el modelo fraccional en lo concerniente a la precisión que se tiene en la caracterización de la respuesta al lograr reproducir los efectos más mínimos oscilatorios.
5Conclusiones
Se comprobó la utilidad del modelo de Zener y por añadidura el del cálculo fraccional en la modelación de tejidos blandos particularmente en un segmento de arteria, los cuales se caracterizaron a través de un modelo constitutivo que obedece las leyes del material viscoelástico.
Se determinó el orden de la derivada fraccionaria que mejor ajusta el modelo aquí aplicado al conjunto de datos, obteniendo un valor que si bien se aproxima a las investigaciones realizadas actualmente por Craiem [18] para las arterias, este es aún más cercano al reportado por Magin y Nagehan [21,26] para otros tipos de tejidos blandos. Se comprobó la importancia de la precisión en la determinación del orden fraccional debido a las diferencias que se presentan entre ambos valores de α a pesar de la cercanía de los valores del orden fraccional.
Se aplicaron una serie de pulsos de entrada al modelo, obteniéndose respuestas que dejan ver una dinámica compleja del sistema de ecuaciones diferenciales y una sensibilidad en los cambios en las condiciones iniciales. Esto lleva a pensar en un análisis propio de la dinámica del sistema y explorar en esta dimensión fractal la viabilidad del caos, como se hace referencia en los artículos de Magin y Awrejcewicz [26,27]. Si bien el modelo de Zener presenta ventajas sobre el de orden entero, este aún resulta limitado ya que en el mismo no se consideran los efectos metabólicos presentes en el cuerpo humano, lo que genera una reconstitución y modificación mecánica del material. De igual manera es necesario escalar este sistema unidimensional a un sistema tridimensional que permita caracterizar los efectos de anisotropía presente en la mayor parte de los materiales biológicos, así como la posibilidad de generar una subrutina de elementos finitos que incluya la formulación fraccional como la desarrollada por Agrawal [28] y no solamente se incluya a partir de la obtención de parámetros del dominio de la frecuencia como en el caso de esta investigación.
Los autores agradecen el apoyo otorgado por el Tecnológico Nacional de México, el Instituto Tecnológico Superior de Cajeme. Al grupo de biomecánica perteneciente al Instituto Superior Politécnico José Antonio Echeverría de la Habana, Cuba.