Вязкоупругая модель склеральной оболочки глаза

МАГИСТЕРСКАЯ РАБОТА
Автор работы: К.П. Фролова
Научный руководитель: к.ф-м.н., доцент Е.Н. Вильчевская

Введение и мотивация работы

В первом приближении глаз имеет шаровидную форму. Он заполнен прозрачной гелеобразной субстанцией, называемой стекловидным телом. Во внешней фиброзной оболочке глаза можно выделить наружную фиброзную, среднюю сосудистую и внутреннюю сетчатую оболочки. Наружная фиброзная (корнеосклеральная) оболочка выполняет защитную функцию и обуславливает форму глаза. Она состоит из передней прозрачной части - роговицы, - и задней непрозрачной части – склеры, обладающих различными радиусами кривизны и биомеханическими характеристиками. Склера занимает 93% внешней фиброзной оболочки глаза человека, поэтому в задачах, связанных с определением формы и изменением объема глазного яблока под действием внутриглазного давления (ВГД), биомеханические свойства склеры играют решающую роль, и поведение глаза допустимо моделировать поведением его склеральной оболочки.

На экспериментальных кривых, характеризующих изменение ВГД в течение нескольких минут после введения интравитреальной инъекции - инъекции внутрь стекловидного тела, - наблюдается его резкий скачок непосредственно после инъекции, вызванный увеличением объема, а затем спад до некоторого постоянного значения. В большинстве существующих источников литературы данный спад объясняется наличием интенсивного оттока внутриглазной жидкости из нагруженного глаза. В связи с уменьшением объема происходит уменьшение ВГД. Однако известно, что склере присуща вязкоупругая реакция на приложенную нагрузку. Непосредственное измерение вязкости склеры вызывает технические сложности и, в связи с этим, в литературе отсутствуют сведения о соответствующих параметрах вязкости. Это приводит к тому, что в большинстве существующих моделей вязкие свойства склеры игнорируются, а поведение склеры при нагрузках предполагается чисто упругим. В данной работе релаксация напряжений (спад ВГД) в глазу после введения интравитреальной инъекции объясняется не только наличием оттока внутриглазной жидкости из нагруженного глаза, но и вязкими свойствами склеральной оболочки глаза. В рамках данной работы будет предложен метод определения коэффициента сдвиговой вязкости склеры, заключающийся сравнении результатов математического моделирования и указанных в литературе экспериментальных данных, основанных на дискретном измерении ВГД в течение нескольких минут после интравитреальной инъекции. Будут рассмотрены различные варианты постановки граничных условий. В первом случае будет предполагаться, что введенный при инъекции дополнительный объем жидкости сохраняется в стекловидном теле на протяжении времени проведения эксперимента, релаксация напряжений будет объясняться только наличием вязких свойств склеры. Во втором случае будет учитываться отток внутриглазной жидкости, релаксация напряжений будет объясняться наличием обоих факторов: и наличием вязких свойств склеры, и наличием интенсивного оттока внутриглазной жидкости из нагруженного глаза. Будет определяться значение коэффициента сдвиговой вязкости, при котором отклонение теоретических данных от экспериментальных минимально.

Цели исследования

• Смоделировать вязкоупругое поведение склеры после интравитреальной инъекции — инъекции внутрь стекловидного тела;
• Определить значение коэффициента сдвиговой вязкости склеры путем сравнения результатов математического моделирования и указанных в литературе экспериментальных данных, основанных на дискретном измерении внутриглазного давления (ВГД) в течение нескольких минут после интравитреальной инъекции.

Постановка задачи

Моделируемый в рамках данной работы эксперимент основан на дискретном измерении ВГД в течение нескольких минут после введения интравитреальной инъекции в объеме 0,05 мл. Точки на экспериментальной кривой соответствуют средним значениям ВГД для 34 пациентов. Для контроля состояния пациентов ВГД измерялось также и на парном не вакцинированном глазу. Задача моделируется вязкоупругим сферическим слоем с внутренним радиусом [math]R_1[/math] и внешним радиусом [math]R_2[/math] при центральносимметричной нагрузке: внешнее давление отсутствует, на внутреннем радиусе заданы перемещения, учитывающие величину дополнительного объема жидкости, введенного при инъекции. Задание нулевого давления на внешнем радиусе объясняется тем, что по определению ВГД есть разница между атмосферным давлением и давлением в глазу. Материал склеры предполагается линейным трансверсально – изотропным. Задача решается в рамках трехмерной теории линейной вязкоупругости.

Материалы и методы

Математическая модель
Система уравнений
Рассматриваемая задача является квазистатической, следовательно, уравнение движения сферического слоя сводится к уравнению равновесия:
[math]\nabla \cdot {\pmb{\sigma }} = 0,[/math]
где [math]\pmb{\sigma }[/math] - тензор напряжений. В силу симметрии задачи данное уравнение в координатном виде в сферической системе координат определяется следующим выражением:
[math]\dfrac{{\partial {\sigma _{rr}}}}{{\partial r}} + 2\dfrac{{{\sigma _{rr}} - {\sigma _{\varphi \varphi }}}}{r} = 0.[/math]
Для получения определяющих соотношений воспользуемся реологической моделью Кельвина – Фойгта, достаточно хорошо описывающей поведение вязкоупругих твердых тел. Данная модель предполагает суммирование упругих и вязких напряжений и равенство упругих и вязких деформаций в теле. При этом для получения единственного решения задачи нам необходимо включать в уравнения лишь один неизвестный параметр вязкости. Итак, определяющие соотношения:
[math]{\pmb{\sigma }} ={{}^4}{\pmb{C}}:{\pmb{\varepsilon }}+ 2\eta {\pmb{\dot e}},[/math]
где [math]{{}^4}{\pmb{C}}[/math] - тензор жесткости четвертого ранга, [math]\eta[/math] - коэффициент сдвиговой вязкости, [math]\pmb{e}[/math] - девиавтор тензора деформаций.
Ненулевые компоненты тензора деформаций связаны с перемещениями следующим образом:
[math]{\varepsilon _{rr}} = \dfrac{{d{u_r}}}{{dr}}\;,\;\; {\varepsilon _{\varphi \varphi }} = {\varepsilon _{\theta \theta }} = \dfrac{{{u_r}}}{r}.[/math]
Безразмерная постановка задачи
Задача решается в безразмерной постановке. Введем следующие безразмерные параметры:
безразмерную координату по радиусу [math]\beta = \dfrac{{{R_1}}}{{{R_2}}} \le x = \dfrac{r}{{{R_2}}} \le 1[/math], безразмерное перемещение [math]u = \dfrac{{{u_r}}}{{{R_2}}}[/math], безразмерную компоненту тензора напряжений [math]{\sigma _{ij\dim }} = \dfrac{{{\sigma _{ij}}}}{{{E_{\theta \theta }}}}[/math], безразмерный модуль Юнга [math]\xi = \dfrac{{{E_{rr}}}}{{{E_{\theta \theta }}}}[/math] и безразмерное время [math]\tau = \dfrac{{{E_{\theta \theta }}}}{\eta }t.[/math]
Тогда уравнение равновесия в перемещениях в безразмерном виде примет следующий вид:
[math]\dfrac{{{\xi ^2}(1 - {\nu _{\theta \varphi }})}}{{\xi (1 - {\nu _{\theta \varphi }}) - 2\nu _{r\theta }^2}}\left[ {\dfrac{{{\partial ^2}u}}{{\partial {x^2}}} + \dfrac{2}{x}\dfrac{{\partial u}}{{\partial x}} - \dfrac{{2(1 - {\nu _{r\theta }})}}{{\xi (1 - {\nu _{\theta \varphi }})}}\dfrac{u}{{{x^2}}}} \right] + \dfrac{4}{3}\left[ {\dfrac{{{\partial ^2}\dot u}}{{\partial {x^2}}} + \dfrac{2}{x}\dfrac{{\partial \dot u}}{{\partial x}} - 2\dfrac{{\dot u}}{{{x^2}}}} \right] = 0.[/math]
ВГД может быть определено как радиальное напряжение на внутренней границе сферического слоя, взятое с обратным знаком, поскольку внешняя нормаль к внутренней границе тела направлена вовнутрь, а ВГД стремится увеличить объем тела. Итак, безразмерное ВГД определяется следующим выражением: [math]IO{P_{\dim}}(\tau) = {\left. { - {\sigma _{xx}}}(\tau) \right|_{x = \beta }}.[/math]
Безразмерное радиальное напряжение как функция перемещений задается следующим образом:
[math]{\sigma _{xx}} (\tau) = \dfrac{{2{\nu _{r\theta }}\xi }}{{\xi (1 - {\nu _{\theta \varphi }}) - 2\nu _{r\theta }^2}}\dfrac{u}{x} + \dfrac{{(1 - {\nu _{\theta \varphi }}){\xi ^2}}}{{\xi (1 - {\nu _{\theta \varphi }}) - 2\nu _{r\theta }^2}}\dfrac{{\partial u}}{{\partial x}} + \dfrac{4}{3}\left[ {\dfrac{{\partial{\dot u}}}{{\partial x}} - \dfrac{{\dot u}}{x}} \right].[/math]
Граничные условия
В рамках данной работы используются два типа граничного условия на внутреннем радиусе. В первом случае предполагается, что объем глазного яблока, включающий в себя дополнительный объем жидкости [math]{\Delta V}[/math], введенный при инъекции внутрь стекловидного тела, и объем глазного яблока [math]V_0^{eye}[/math] до нагружения, сохраняется на протяжении времени проведения эксперимента. При этом напоминаем, что моделирование глазного яблока сводится к моделированию склеральной оболочки глаза. Граничные условия задаются следующим образом:
[math]{\sigma _{xx}}(\tau )\left| {_{x = 1}} \right. = 0\;,\;\; u(\tau )\left| {_{x = \beta }} \right. = {u_0}H(\tau ),[/math]
где [math]{u_0} \approx \dfrac{{\Delta V}}{{4\pi {R_2}R_1^2}}[/math], [math]H(\tau)[/math] - единичная степенная функция Хевисайда.
При задании второго типа граничного условия на внутреннем радиусе учитывается гидродинамика внутриглазной жидкости, в частности, более интенсивный отток внутриглазной жидкости из нагруженного глаза. Для определения величины изменения объема глазного яблока, обусловленного притоком и оттоком внутриглазной жидкости, обратимся к методу тонографии. Данный метод оперирует скоростью изменения объема глазного яблока. Текущий объем глаза определяется следующим образом:
[math]V(t) = V(0)+\int\limits_{\tilde t = 0}^{\tilde t = t} {\dot V(\tilde t)d\tilde t} = V_0^{eye} + \Delta V + \int\limits_{\tilde t = 0}^{\tilde t = t} {\dot V(\tilde t)d\tilde t}.[/math]
Изменение объема глаза, вызванное гидродинамикой внутриглазной жидкости: [math]V(t) - V(0) = V(t) - V_0^{eye} - \Delta V = 4\pi R_{inj}^2{u_r}(t),[/math] где [math]{R_{inj}} = \sqrt[3]{{R_1^3 + \dfrac{{3\Delta V}}{{4\pi }}}}.[/math] Граничные условия задаются следующим образом:
[math]{\sigma _{xx}}\left( t \right)\left| {_{x = 1}} \right. = 0,u\left( t \right)\left| {_{x = \beta }} \right. = {u_0} + \dfrac{{\int\limits_{\tilde t = 0}^{\tilde t = t} {\dot V\left( {\tilde t} \right)d\tilde t} }}{{4\pi {R_2}R_{inj}^2}}.[/math]
Дальнейшее использование безразмерного времени затруднительно, поскольку в граничном условии на внутреннем радиусе присутствует интеграл по времени. Для определения величины [math]\int_{\tilde t = 0}^{\tilde t = t} {\dot V(\tilde t)d\tilde t}[/math] воспользуемся данными тонографии - метода измерения и регистрации ВГД, позволяющего определить интенсивность оттока внутриглазной жидкости. В рамках данного метода предполагается, что скорость изменения объема глазного яблока зависит от скоростей притока и оттока внутриглазной жидкости. При этом получается, что интересующий нас интеграл определяется следующим образом:
[math]\int\limits_{\tilde t = 0}^{\tilde t = t} {\dot V\left( {\tilde t} \right)d\tilde t} = \int\limits_{\tilde t = 0}^{\tilde t = t} {\left( {F\left( {\tilde t} \right) - R\left( {\tilde t} \right)} \right)d\tilde t},[/math]
где [math]F, R[/math] - скорости притока и оттока внутриглазной жидкости соответственно.
Величина оттока определяется следующим гидравлическим соотношением: [math]R = C\left( {P - {P_e}} \right),[/math] где [math]C[/math] - коэффициент легкости оттока внутриглазной жидкости, [math]P[/math] - аппроксимирующая функция ВГД, [math]{P_e}[/math] - давление в эписклеральных венах. Основная задача тонографии состоит в том, чтобы оценить величину коэффициента легкости оттока. При стандартной обработке данных тонографии данная величина оценивается на основе представления тонограммы линейной функцией, причем для обработки используются только начальная и конечная точки тонограммы. В работе Г. Любимова предложен модифицированный алгоритм обработки данных тонографического исследования, предусматривающий оценку тонографической кривой, основанную на использовании не только начального и конечного значений ВГД, но и его промежуточных значений. При этом рассматривается несколько схем обработки данных тонограммы. При стандартной схеме делаются следующие предположения: [math]{P_e} - {P_{e0}} = 1.25[/math], [math]C = {C_0} = const[/math], [math]F = {F_0} = const, [/math] где индекс "0" соответствует параметрам в ненагруженном состоянии.
При 1-ой модифицированной схеме делаются следующие предположения: [math]{P_e} - {P_{e0}} = 1.25[/math], при этом [math]F = {R_{st}} \ne {F_0}[/math]. При 2-ой модифицированной схеме предполагается, что [math]F = {R_{st}} = {F_0}[/math], при этом [math]{P_e} - {P_{e0}} \ne 1.25[/math].
Метод преобразования Лапласа
Для получения решения для перемещений и радиальных напряжений воспользуемся методом преобразования Лапласа, связывающим функцию [math]f\left( {x,\tau } \right)[/math] действительного переменного (оригинал), с функцией [math]\bar f\left( {x,s} \right)[/math] комплексного переменного (изображением). Изображение функции - оригинала по Лапласу определяется следующим выражением:
[math]\bar f\left( {x,s} \right) = \exp \left( { - s\tau } \right)f\left( {x,\tau } \right)d\tau,[/math]
где [math]s[/math] - комплексная переменная, [math]s = c + i\omega ,c \gt {s_0}[/math].
Преимуществом использования метода Лапласа в рамках данной работы является замена оператора дифференцирования по времени оператором умножения на комплексную переменную:
[math]\dot f(x,\tau ) \to \int\limits_0^\infty {{e^{ - s\tau }}} \dot f(x,\tau )d\tau = s\bar f(x,s) - f(x,0).[/math]
Постановка начального условия задачи затруднительна в силу сложности экспериментального определения перемещений в склере в начальный момент времени. В рамках данной работы мы будем предполагать, что в начальный момент времени перемещения зависят от радиальной координаты линейно. Тогда слагаемые в преобразованных по методу Лапласа выражениях для уравнения равновесия и радиальных напряжений, содержащие величину перемещений в начальный момент времени, сократятся. Итак, в пространстве Лапласа оперируем следующими дифференциальными уравнениями, содержащими только производные по координате:
[math]\left[ {\dfrac{{{\xi ^2}(1 - {\nu _{\theta \varphi }})}}{{\xi (1 - {\nu _{\theta \varphi }}) - 2\nu _{r\theta }^2}} + \dfrac{4}{3}s} \right]\left( {\dfrac{{{\partial ^2}\bar u}}{{\partial {x^2}}} + \dfrac{2}{x}\dfrac{{\partial \bar u}}{{\partial x}} - 2\dfrac{{\bar u}}{{{x^2}}}} \right) - \dfrac{{2\xi [1 - {\nu _{r\theta }} - \xi (1 - {\nu _{\theta \varphi }})]}}{{\xi (1 - {\nu _{\theta \varphi }}) - 2\nu _{r\theta }^2}}\dfrac{{\bar u}}{{{x^2}}} = 0;[/math]
[math]\overline {{\sigma _{xx}}} = \left[ {\dfrac{{2{\nu _{r\theta }}\xi }}{{\xi (1 - {\nu _{\theta \varphi }}) - 2\nu _{r\theta }^2}} - \dfrac{4}{3}s} \right]\dfrac{{\bar u}}{x} + \left[ {\dfrac{{(1 - {\nu _{\theta \varphi }}){\xi ^2}}}{{\xi (1 - {\nu _{\theta \varphi }}) - 2\nu _{r\theta }^2}} + \dfrac{4}{3}s} \right]\dfrac{{d\bar u}}{{dx}}.[/math]
Для решения данных уравнений необходимо воспользоваться граничными условиями, преобразованными по правилам Лапласа. Решение для изображений функций перемещений и радиальных напряжений ищется в пакете символьной математики WolframMathematica с помощью команд Solve и DSolve.
Численные алгоритмы обратного преобразования Лапласа
Для определения оригиналов функций перемещений и радиальных напряжений необходимо воспользоваться формулой обращения (интегралом Бромвича) – формулой обратного преобразования Лапласа. Она имеет следующий вид:
[math]f(x,\tau ) = \dfrac{1}{{2\pi i}}\int\limits_C {\bar f(x,s)} \exp (s\tau )ds,\tau \gt 0,[/math]
где [math]C[/math] - контур от [math]c - i\infty[/math] до [math]c + i\infty[/math]
Применим данное обращение к найденным изображениям функций перемещений и напряжений. Для трансверсально – изотропного случая найденные функции имеют сложную зависимость от набора параметров и, в связи с этим, получение аналитического решения на основе интеграла Бромовича не представляется возможным. Обратимся к численным методам обратного преобразования Лапласа. При решении задачи в рамках данного численного подхода применяется квадратурная формула типа Гаусса:
[math]f(x,\tau) \approx {f_n}(x,\tau ) \equiv {f_{n,a,K}}(x,\tau ) = \frac{1}{\tau }\sum\limits_{j = 1}^n {{K_j}\bar f} (x,\frac{{{a_j}}}{\tau }),[/math]
где [math]\bf {a},\,\bf {K}[/math] - векторы, называемые узлами и весами соответственно.
Данное выражение можно получить из интеграла Бромовича путем разложения экспоненциальной функции в ряд Маклорена и дальнейшего применения аппроксимации Паде к полученной таким образом степенной функции. Существует несколько алгоритмов определения весов и узлов квадратурной формулы. Одним из наиболее известных алгоритмов является алгоритма Закина, используемая при этом квадратурная формула выражается следующим образом: [math]f_{n,a,K}^{Zakian}\left( {x,\tau } \right) = 2{f_{n,a,K}}\left( {x,\tau } \right)[/math]. В качестве альтернативных алгоритмов мы рассмотрим алгоритмы Джеффресона и Чоу для 10-ти и 15-ти членов разложения ряда в квадратурной формуле, задающей преобразование Лапласа.

Результаты исследования

Мы получили реальные функции перемещений и радиальных напряжений, зависящие от времени, путем применения численного обратного преобразования Лапласа к полученным нами ранее функциям изображений для перемещений и радиальных напряжений. Мы исследовали результаты, полученные на основе применения трех перечисленных в данной работе численных алгоритмов, а именно: алгоритма Закиана, алгоритма Джеффресона и Чоу для 10-ти и 15-ти членов разложения ряда в квадратурной формуле, задающей численное обратное преобразование Лапласа. Коэффициент сдвиговой вязкости склеры был определен при двух типах граничного условия на внутреннем радиусе: в первом случае предполагалось, что объем глаза не меняется в течение времени проведения эксперимента, во втором учитывался отток внутриглазной жидкости. Зависимость скорости изменения объема внутриглазной жидкости от времени определялась на основании нескольких схем обработки данных тонографического исследования, предложенных в работе Г. Любимова. При численном решении задачи мы использовали следующие значения параметров: [math]{{R_1}}=11.75[/math] мм, [math]{{R_2}}=12.25[/math] мм, [math]{{E_{22}}}=14.3[/math] МПа, [math]{{E_{11}}}=0.01{{E_{22}}}[/math], [math]{{\nu _{12}}}=0.01[/math], [math]{{\nu _{23}}}=0.45[/math], указанные в работах С.М. Бауэр.
Результаты, полученные в предположении постоянства объема глаза на протяжении времени проведения эксперимента
Обсудим результаты, полученные в предположении, что объем глаза не меняется на протяжении времени проведения эксперимента. При этом будем пользоваться стандартной схемой обработки данных тонографического исследования. В данном случае мы используем все безразмерные величины, перечисленные ранее в настоящей работе. Отметим, что в данном случае безразмерное радиальное напряжение не зависит от коэффициента сдвиговой вязкости в явном виде. Данный коэффициент не входит в уравнение равновесия и в выражение для радиальной компоненты тензора напряжений и не задается в качестве известного параметра задачи. Безразмерное радиальное напряжение является функцией безразмерного времени, которое, в свою очередь, связано с коэффициентом сдвиговой вязкости через размерное время. Так, коэффициент сдвиговой вязкости может быть определен по следующей формуле:
[math]\eta = \dfrac{{{E_{\theta \theta }}t}}{\tau }.[/math]
Рассмотрим экспериментальные данные, полученные независимо в работах К. Котляра и Б. Першина. Экспериментальные данные, полученные в первой работе: [math]ВГД = 6266[/math] Па [math](ВГ{Д_{безр}} \cdot 10^{-6} = 438.2), t = 10[/math] с; [math]ВГД = 4533[/math] Па [math](ВГ{Д_{безр}} \cdot 10^{-6} = 317.0), t = 120[/math] с; [math]ВГД = 3466[/math] Па [math](ВГ{Д_{безр}} \cdot 10^{-6} = 242.4), t=300[/math] с; [math]ВГД =2800[/math] Па [math](ВГ{Д_{безр}} \cdot 10^{-6} = 195.8), t = 500[/math] с. Экспериментальные данные, полученные во второй работе: [math]ВГД = 8733[/math] Па [math](ВГ{Д_{безр}} \cdot 10^{-6} = 610.7), t = 0[/math] с; [math]ВГД = 6853[/math] Па [math](ВГ{Д_{безр}} \cdot 10^{-6} = 479.2), t = 60[/math] с; [math]ВГД = 4800[/math] Па [math](ВГ{Д_{безр}} \cdot 10^{-6} = 335.6), t=180[/math] с; [math]ВГД =2800[/math] Па [math](ВГ{Д_{безр}} \cdot 10^{-6} = 195.8), t = 300[/math] с. Определяем, при каком значении координаты по безразмерному времени полученная нами функция безразмерного радиального напряжений равна по модулю перечисленным экспериментальным значениям ВГД. Оказывается, что диапазон значений ВГД, получаемый теоретически, значительно меньше диапазона экспериментальных значений. Так, диапазон, полученный теоретически: [math]IO{P_{\dim}}\cdot 10^{-6}\in 285-430[/math], - тогда, как диапазон экспериментальных значений: [math]IO{P_{\dim}}\cdot 10^{-6}\in 195.8-438.2[/math] для данных, представленных в работе К. Котляра и [math]IO{P_{\dim}}\cdot 10^{-6}\in 195.8-610.7[/math] для данных, представленных в работе Б. Першина. Таким образом, мы можем использовать только данные для второй точки эксперимента, рассмотренного в работе К. Котляра (ВГД=4533 Па, t=120 с), и данные для третьей точки эксперимента, рассмотренного в работе Б. Першина (ВГД=4800 Па, t=180 с), для определения координаты по безразмерному времени. Значения коэффициента сдвиговой вязкости склеры для двух указанных экспериментальных точек при применении численного алгоритма Закиана и численных алгоритмов Джеффресона и Чоу для 10-ти и 15-ти членов разложения ряда в квадратурной формуле, задающей численное обратное преобразование Лапласа, представлены в таблице:

Экспериментальная точка	Алгоритм Закиана	Алгоритм Д&Ч, n=10	Алгоритм Д&Ч, n=15
ВГД=4533 Па, t=120 с	[math]\eta = 12.8[/math] МПа [math]\cdot[/math] с	[math]\eta = 12.8[/math] МПа [math]\cdot[/math] с	[math]\eta = 12.9[/math] МПа [math]\cdot[/math] с
ВГД=4800 Па, t=180 с	[math]\eta = 27.8[/math] МПа [math]\cdot[/math] с	[math]\eta = 27.7[/math] МПа [math]\cdot[/math] с	[math]\eta = 27.8[/math] МПа [math]\cdot[/math] с

Кривые релаксации напряжений, построенные при использовании трех различных численных алгоритмов, совпадают при указанном на рисунке справа масштабе. При построении графиков мы использовали метод Закиана и соответствующие ему значения коэффициента сдвиговой вязкости. Основной вывод, который мы можем сделать, анализируя график, изображенный на рисунке, заключается в том, что при использовании предположения о постоянстве объема глаза на протяжении времени проведения эксперимента построенная нами теория не достаточно хорошо описывает экспериментальные данные как количественно, так и качественно. Теоретические кривые обладают менее резким спадом, чем экспериментальные кривые. Так, релаксация напряжений не может быть объяснена только наличием вязкости.


Результаты, полученные в предположении наличия интенсивного оттока внутриглазной жидкости из нагруженного глаза
Перейдем к обсуждению результатов, полученных при учете гидродинамики внутриглазной жидкости. В данном случае мы используем размерное время (в секундах) при построении уравнения равновесия и выражения для радиальной компоненты тензора напряжений в связи с тем, что в граничном условии для перемещений присутствует интеграл по времени. Коэффициент сдвиговой вязкости склеры в данном случае входит в решаемые уравнения в явном виде, а его величину необходимо задавать в качестве известного параметра при численных вычислениях. Для определения оптимального значения коэффициента сдвиговой вязкости воспользуемся методом бисекции, основанном на первой теореме Больцано – Коши: [math]\Phi (\eta ) \in C[{\eta _a},{\eta _b}],\Phi ({\eta _a}) \cdot \Phi ({\eta _b}) \lt 0 \Rightarrow \exists {\eta _c} \in [{\eta _a},{\eta _b}]:\Phi ({\eta _c}) = 0[/math].В качестве функции [math]\Phi (\eta )[/math] будем рассматривать:
[math]\Phi \left( \eta \right) = \left( {IO{P_{\exp eriment(\dim )}} - \left( { - {{\left. {{\sigma _{xx}}} \right|}_{x = \beta }}} \right)} \right)\left| {_{t = {t_{{\rm{experiment}}}}}} \right..[/math]
Метод бисекции позволяет нам найти корни данного уравнения для всех экспериментальных точек, представленных в работе К. Котляра, кроме первой точки, характеризующей скачок ВГД непосредственно после инъекции. При рассмотрении экспериментальных значений ВГД, представленных в работе Б. Першина оказалось, что корней уравнения не существует. Значения коэффициента сдвиговой вязкости склеры, полученные итерационным методом при использовании численного алгоритма Закиана и численных алгоритмов Джеффресона и Чоу для 10-ти и 15-ти членов разложения ряда в квадратурной формуле, задающей численное обратное преобразование Лапласа, представлены в таблице ниже. Значения определены с точностью до второго знака после запятой.

Экспериментальная точка	Алгоритм Закиана	Алгоритм Д&Ч, n=10	Алгоритм Д&Ч, n=15
ВГД=4533 Па, t=120 с	[math]\eta = 57.1[/math] МПа [math]\cdot[/math] с	[math]\eta = 57.1[/math] МПа [math]\cdot[/math] с	[math]\eta = 57.1[/math] МПа [math]\cdot[/math] с
ВГД=3466 Па, t=300 с	[math]\eta = 58.6[/math] МПа [math]\cdot[/math] с	[math]\eta = 59.1[/math] МПа [math]\cdot[/math] с	[math]\eta = 58.3[/math] МПа [math]\cdot[/math] с
ВГД=2800 Па, t=500 с	[math]\eta = 33.8[/math] МПа [math]\cdot[/math] с	[math]\eta = 35.3[/math] МПа [math]\cdot[/math] с	[math]\eta = 33.0[/math] МПа [math]\cdot[/math] с

Функции зависимости ВГД от времени представлены на рисунке справа. Кривые релаксации напряжений, построенные при использовании трех различных численных методов, совпадают при указанном на рисунке масштабе. При построении графиков мы использовали метод Закиана и соответствующие ему значения коэффициента сдвиговой вязкости. Как видно из рисунка, значение коэффициента сдвиговой вязкости склеры влияет на характер спада кривой, соответствующей релаксации напряжений. При большем коэффициенте сдвиговой вязкости наблюдается более продолжительная релаксация напряжений. В случае, когда коэффициент сдвиговой вязкости равен нулю, что соответствует отсутствию вязкостных свойств склеры в рамках принятой модели, теория не достаточно удовлетворительно описывает эксперимент. Это означает, что релаксацию напряжений после введения интравитреальной инъекции нельзя объяснить только наличием оттока внутриглазной жидкости.
Сравнивая результаты, полученные при двух типах граничного условия для перемещений на внутреннем радиусе, мы можем отметить, что значение коэффициента сдвиговой вязкости склеры меньше в случае, когда объем глаза предполагается фиксированным на протяжении времени проведения эксперимента. Также мы можем убедиться в том, что необходимо принимать во внимание оба фактора: и наличие вязкости склеры, и учет интенсивного оттока внутриглазной жидкости из нагруженного глаза, - для лучшего удовлетворения теорией экспериментальных данных.
Результаты, полученные при разных схемах обработки тонограммы

Обсудим результаты, полученные при использовании модифицированных схем обработки данных тонографического исследования. Как мы уже показали, численные алгоритмы обратного преобразования Лапласа, рассматриваемые в рамках данной работы, дают практически одинаковые результаты. В связи с этим далее мы будем использовать только алгоритм Закиана для получения оригиналов функций перемещений и радиальных напряжений, зависящих от времени. Рассмотрим результаты, полученные в предположении зависимости величины скорости притока внутриглазной жидкости от приложенной нагрузки, в то время как давление в эписклеральных венах предполагается практически не зависимым от нагружения. Значение коэффициента сдвиговой вязкости: [math]\eta = 17.8 МПа\cdot с[/math]. Функция зависимости ВГД от времени представлена на рисунке слева. Видно, что кривые, характеризующие изменение ВГД с течением времени, стремятся к асимптоте, проходящей выше двух последних экспериментальных точек. Следовательно, применение обсуждаемой схемы обработки данных тонографического исследования приводит к результатам, хуже согласующимся с экспериментальными данными, чем результаты, полученные путем применения стандартной схемы.

Рассмотрим результаты, полученные в предположении равенства притоков внутриглазной жидкости в нагруженном и ненагруженном глазу, в то время как давление в эписклеральных венах не фиксируется с помощью заданного соотношения [math]{P_e} - {P_{e0}} = 1.25[/math] мм.рт.ст. Значение коэффициента сдвиговой вязкости, полученное методом бисекции: [math]\eta = 67.6 МПа \cdot с[/math]. Функция зависимости ВГД от времени представлена на рисунке справа. Видно, что полученные теоретически кривые, характеризующие изменение ВГД с течением времени, имеют более резкий спад, чем экспериментальная зависимость. При этом не наблюдается значение ВГД, к которому стремились бы теоретические кривые, что не соответствует действительности, поскольку интенсивный отток должен происходить только из нагруженного глаза, т.е. до того момента, пока дополнительный объем жидкости, введенный при инъекции, не вытечет. После этого ВГД должно прийти в норму, а не свестись к нулю. Следовательно, применение обсуждаемой схемы обработки данных тонографического исследования приводит к результатам, хуже согласующимся с экспериментальными данными, чем результаты, полученные путем применения стандартной схемы.
Результаты, полученные при разных комбинациях значений модулей Юнга

Исследуем влияние модулей упругости материала склеры на значение коэффициента сдвиговой вязкости. Данное исследование актуально в силу того, что в разных источниках литературы приводятся разные значения модулей Юнга в направлении оси симметрии и в плоскости изотропии. Рассмотрим комбинацию модулей Юнга, близкую к рассматриваемой нами в рамках данной работы: [math]{E_{\theta \theta }} = 17.4[/math] МПа, [math]{E_{rr}} = 0.01{E_{\theta \theta }}[/math]. Применяя метод бисекции, можно получить значения коэффициента сдвиговой вязкости склеры для второй экспериментальной точки, приведенной в работе К. Котляра (ВГД=4533 Па, t=120 с) и для третьей экспериментальной точки, приведенной в работе Б. Першина (ВГД=4800 Па, t=180 с): [math]\eta = 17.6[/math] МПа [math]\cdot[/math] с, [math]\eta = 57.2[/math] МПа [math]\cdot[/math] с соответственно. Функция зависимости ВГД от времени представлена на рисунке слева. Можно сделать вывод о том, что коэффициент сдвиговой вязкости склеры зависит от значений модулей Юнга, при этом чем больше значения модулей Юнга, тем меньше значение сдвиговой вязкости.
Функции зависимости ВГД от времени при некоторых других комбинациях модулей Юнга представлены на рисунках ниже.

Рисунок слева построен при следующей комбинации модулей Юнга: [math]{E_{\theta \theta }} = 41.0[/math] МПа, [math]{E_{rr}} = 0.01{E_{\theta \theta }}[/math]. Рисунок справа построен при следующей комбинации модулей Юнга: [math]{E_{\theta \theta }} = 20.0[/math] МПа, [math]{E_{rr}} = 0.5[/math] МПа.
Можно заключить, что использование комбинации модулей Юнга [math]{E_{\theta \theta }} = 14.3[/math] МПа, [math]{E_{rr}} = 0.01{E_{\theta \theta }}[/math], предлагаемой в ряде работ, в. т.ч. в работах С.М. Бауэр, оптимально, т.к. в данном случае теория лучше всего согласуется с экспериментальными данными, приведенными в работе К. Котляра.

Выводы

При учете обоих факторов: наличия вязкости склеральной оболочки глаза и интенсивного оттока внутриглазной жидкости из нагруженного глаза, - теория лучше согласуется с экспериментом, чем при учете одного из них. Таким образом, релаксацию напряжений невозможно описать качественно наличием только лишь вязкоупругого поведения склеры при нагружении или наличием гидродинамики внутриглазной жидкости. При этом оказалось, что значение коэффициента сдвиговой вязкости склеры меньше в случае предположения постоянства объема глаза на протяжении времени проведения эксперимента. Величина коэффициента сдвиговой вязкости влияет на характер релаксации напряжений: чем больше коэффициент сдвиговой вязкости, тем более продолжительна релаксация напряжений. Выбор численного алгоритма обратного преобразования Лапласа при рассмотрении алгоритма Закиана, а также алгоритмов Джеффресона и Чоу для 10-ти и 15-ти членов разложения ряда в квадратурной формуле, задающей численное обратное преобразование Лапласа, не оказывает значительного влияния на результат. При этом стоит отметить, что использование алгоритма Закиана наименее трудоемко с вычислительной точки зрения в силу использования наименьшего числа членов разложения ряда в квадратурной формуле. Традиционные предположения тонографии приводят к результатам, лучше согласующимся с экспериментальными данными, чем предположения модифицированных схем обработки данных тонографического исследования. Такой вывод может объясняться тем, что при попытке описать одни физические механизмы работы гидродинамической системы глаза необходимо параллельно учитывать ряд других процессов и факторов. Так, если принимать, что часть параметров, характеризующих состояние глаза, изменяется при нагружении глаза, вполне вероятно, что нужно учитывать сам характер нагружения и величину прикладываемой нагрузки. В рамках же рассматриваемой модели в работе Г. Любимова принято, что ряд параметров меняется при нагружении, но при этом не зависит явно от величины ВГД и приложенной нагрузки, т.е. меняется скачкообразно и затем сохраняется на протяжении времени проведения эксперимента. При этом предположение о том, что скорость притока внутриглазной жидкости меняется при нагружении, а давление в эписклеральных венах остается практически неизменным, приводит к меньшему значению коэффициента сдвиговой вязкости склеры и наблюдению меньшего спада ВГД с течением времени. Предположение о том, что приток жидкости не зависит от ВГД и величины нагрузки, тогда как давление в эписклеральных венах предполагается заметно отличным от соответствующего значения в ненагруженном глазу, приводит к наблюдению резкого спада ВГД с течением времении, что качественно не согласуется с экспериментальными данными. Найденные значения коэффициента сдвиговой вязкости в этом случае больше. Коэффициент сдвиговой вязкости склеры зависит от значений модулей Юнга в направлении оси симметрии и плоскости изотропии в трансверсально – изотропном материале. При допустимом варьировании значений модулей оказывается, что коэффициент сдвиговой вязкости тем больше, чем меньше модули Юнга. При этом использование комбинации модулей Юнга, предложенной в работах С.М. Бауэр, является наиболее удовлетворительной при сопоставлении теоретических данных с экспериментальными.

Список литературы

• D. Ljubimova. Biomechanics of the Human Eye and Intraocular Pressure Measurements // Ph.D. thesis - Royal Institute of Technology - 2009.
• Л.А. Карамшина. Модели многослойных оболочек в задачах офтальмологии // Дисс. – Санкт-Петербургский государственный университет - 2011.
• С.А. Регирер. Лекции по биологической механике – Изд. Московского университета – 1980 – 144 с.
• К. Котляр, С. Бауэр, Н. Планге. Клинические и биомеханические аспекты изменения внутриглазного давления после интравитреальных инъекций // Российский общенациональный офтальмологический форум - 2013.
• В. Новацкий. Теория упругости // Пер. с польского – М.: Мир - 1975 - 872 с.
• Р. Кристенсен. Введение в механику композитов // Пер. с англ. – М.: Мир - 1982 – 334 с.
• Е.Н. Иомдина, С.М. Бауэр, К.Е. Котляр. Биомеханика глаза: теоретические аспекты и клинические приложения // М.: Реал Тайм – 2015 – 208 с.
• Г.А. Любимов, И.Н. Моисеева, А.А. Штейн, Е.Н. Иомдина, Л.А. Назаренко. Об оценке величины оттока жидкости из глаза с помощью модифицированного метода тонографии // Российский журнал биомеханики – 2012 – Т.16 - №2 – сс. 8–20.
• Л.А. Игумнов, С.Ю. Литвинчук, А.А. Белов. Численное обращение преобразования Лапласа: Учебно – методическое пособие // Изд. Нижегородского университета – 2010 – 34 с.
• G. Doetsch. Introduction to the Theory and Application of the Laplace Transformation // Springer – 1974 - 327 p.
• J. Abate, W. Whitt. A Unified Framework for Numerically Inverting Laplace Transforms// INFORMS Journal on Computing – 2006 - V.18 - No.4 - pp.408-421.
• G. A. Baker, P. Graves-Morris. Pad´e Approximants, 2nd ed. // Cambridge University Press, Encyclopedia of Mathematics and Its Applications – 1996 - Vol.59 – 764 p.
• V. Zakian. Numerical inversion of Laplace transform // Elec. Lett. – 1969 - No.5 - pp.120-121.
• V. Zakian. Optimisation of numerical inversion of Laplace transforms // Elec. Lett. – 1970 - No.6 - pp.677-679.
• C.P. Jeffreson, E.-P. Chow. Least squares coefficients for a quadrature formula for Laplace transform inversion // Journal of Computational and Applied Mathematics – 1978 – V.4 - №1 - pp. 53-58.
• G. Miller. Least squares approximation of functions using exponentials // Ph.D. thesis - John Hopkins Univ. - 1969.
• С.М. Бауэр, Л.А. Замураев, К.Е. Котляр. Модель трансверсально-изотропного сферического слоя для расчета изменения внутриглазного давления при интрасклеральных инъекциях // Российский журнал биомеханики – 2010 – Т.10 - №2 – сс.43–49.
• Б.С. Першин. Гидродинамический баланс глазного яблока при интравитреальном введении дополнительного объема жидкости // Дисс. – ФГБУ Научно-исследовательский институт глазных болезней РАМН - 2012.
• C.H. Morales. A Bolzano’s theorem in the new millennium // Nonlinear Analysis – 2002 - No. 51 - pp. 679-691.
• Е.Н. Иомдина. Механические свойства тканей человека // Изд. МГУ: Современные проблемы биомеханики – 2006 – Вып.11 – сс.183-200.