Математическое моделирование повреждения глазного яблока при контакте с летящим осколком

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

Причины, вызывающие повреждения глаза, могут иметь механическую, тепловую, химическую и т.п. природу. Большая доля глазных травм (до 70%) происходит на производстве [1]. Это указывает не только на уязвимость глаза, но и на недостаточные меры безопасности, применяемые против таких травм.

Среди производственных глазных травм наиболее опасным являются травмы с внутриглазным проникновением инородного тела. Такие повреждения глазного яблока обычно имеют летальный исход.

Многочисленные клинические данные о последствиях глазных травм, размеры извлекаемых осколков, а также анализ причин и условия возникновения летящих осколков можно считать результатами натурных экспериментов по проникновению инородного тела в глазное яблоко. На основе этих данных с использованием современных методов математического моделирования и компьютерной техники, появляется возможность изучения не только физических и механических условий глазных травм, но и оценки влияния этих травм на последующее функционирование поврежденного глаза.

В приведенном ниже материале моделируются процессы деформирования и разрыва тканей роговой оболочки и склеры глазного яблока человека при контакте с летящим инородным телом. Моделирование процесса проникновения осколка ограничивается разрывом роговицы (или склеры) и началом вхождением осколка в жидкостную среду первичной камеры.

В основе математической постановки задачи лежит подробный учет реального строения глазного яблока и свойств заполняющих его биологических веществ. Численное решение соответствующей краевой задачи строится с использованием комбинированного метода конечных элементов и конечных разностей соответственно для описания геометрии глаза и решения по пространству и на временном слое.

Постановка задачи. Математические модели

Летящее инородное тело (осколок) будем рассматривать в виде абсолютно твердой частицы клиновидной формы. Скорость полета осколка может определяться, например, как произведение угловой скорости вращения ротора инструмента на его радиус, или скоростью воздушного потока, переносящего абразивные частицы.

Глазное яблоко рассматривается в виде сферических оболочек для роговицы и склеры по границе, заполненных жидкообразной средой с внутриглазным давлением 20 мм рт. ст, снаружи она нагружена атмосферным давлением, действующим на роговицу и склеру.

Отметим, что в такой объемной постановке удается учесть все основные составные части глазного яблока и сложную трехмерную картину распространения в нем ударных волн и процессов повреждения.

Для пространственной аппроксимации передней камеры, хрусталика и стекловидного тела были применены пространственные восьмиузловые, а для роговицы и склеры четырехузловые безмоментные конечные элементы.

Представленное таким образом глазное яблоко рассматривается далее как нелинейно деформируемая среда объема V , занимающая в R ³ произвольную, в общем случае многосвязную область D , ограниченную поверхностью S . В этом случае область D включает в себя поверхности раздела в виде соединения разнородных материалов (сред) и возможных дефектов в виде отслоений и повреждений.

Начальную конфигурацию к ₀ области D отнесем к прямоугольной декартовой системе координат x i , i = 1,2,3 (эйлеровы координаты). Движение (деформацию) тела D относительно начальной конфигурации к ₀( D ) в любой произвольный момент времени t >  0 можно описать следующим отображением:

X , = Z ( X k , t ),

X = (X 1,X2,X3) е к0(D), t > 0, к0 : D ^ R3, где Xk - лагранжевы координаты, а Z - непрерывно обращаемое отображение из к0 в R3.

Мерой деформации тела D служит градиент деформации F = d x /d X , который согласно известной теореме о полярном разложении в теории конечных деформаций может быть представлен в виде

F = RU, где R - положительно определенный ортогональный тензор вращения, U -симметричный тензор растяжения.

Градиент скорости частицы тела D равен:

5 v 5 v 5 X t_T i L = — =--= F F ,

5x 5X 5x этот тензор может быть представлен как сумма тензоров скоростей вращения W и растяжения Y:

L = W + Y, где

1 T

Y = 2( L + L T ),

W = 2 ( L¹ - L ) .                                 (3)

Скорость тензора деформации Грина-Сен Венана E = 0,5(U ² - 1) определяется из соотношения:

E = F^TYF .                                 (4)

Среда движется (деформируется) под действием импульсных ударных воздействий, в общем случае, задаваемых или в виде поверхностных сил F_j ( xⁱ , t ) , x i е ^ S_F , t e D_t = (0, т ), т - время действия импульса, или начального распределения скоростей V j ( x i , 0) на части поверхности S V , причем S = S_F и S V и S_F n S V = 0 . Здесь и далее в отношении индексов и операций над индексами используются соглашения, принятые в тензорном исчислении.

Система уравнений импульсной динамики деформируемой среды, описывающая ее механическое поведение, может быть представлена в следующем виде [2]:

р det F = ро, div с + р f = рvt,     на D х Dt                         (5)

e = VSj а j- (p + q V, где р - плотность рассматриваемой среды, верхний знак t означает производную по времени, f - массовые силы, vi = ut - вектор скорости, с, а - тензоры напряжений и деформаций, соответственно, Sy =су - p6ij - девиатор тензора напряжений, p = сii /3 – давление, e – полная энергия на единицу объема среды, V – относительный объем текущего состояния, q – объемная вязкость.

В начальный момент времени t = 0

р = р о , p = p 0 , e = e 0 , x ( X^k ,0) = X^k , x t ( X^k ,0) = v ( X^k ,0),            (6)

где v_i – начальное распределение скоростей.

Приведем граничные условия для деформируемой среды D . Запишем их для известной формы начальной границы S ₀ в момент времени t = 1 ₀, S ₀ = S ( 1 ₀):

x ( X^k , t ) = r i (t ), X^k e S ”,                                  (7)

a j n, = F , (X k , t ), X k e S F ,

(7а)

При наличии внутренней границы контактного разрыва S C e S ⁰ к приведенным условиям (7) добавляются условия контактной непрерывности:

(a. - a2) n. = 0, Xk = X2, X1 e SC x D, где цифрами 1, 2 обозначены тела (или разнородные материалы), находящиеся по обе 0

стороны от границы S_c .

Если при этом возникает возможность относительного движения указанных тел с границами S 0 и S_c ^o2 , первоначально совпадающими с S C , условию отрыва общих точек поверхности в зоне контакта соответствуют нулевые усилия в них, ^{( a}n- ^a 2 jnj = ⁰⁾.

Условию скольжения с трением соответствует ( -u F n/² <  F t ^1,2<ц ^ 1’²), где индексами t и n отмечены касательное и нормальное направления к поверхностям в точке контакта, ц - коэффициент трения, F n ^1,2 = a . n i n . , F¹ "² = a . n i t . , t . - единичный вектор касательной.

Для уравнения состояния деформируемой среды мы ограничимся так называемым квазиакустическим приближением, когда давление является только функцией объема, а ударная адиабата совпадает с изоэнтропой расширения [2]:

p = K | -1- - 1 | , ( V )

где K – объемный модуль упругости

Для сглаживания разрывов в профилях изменения плотности, напряжений и других полевых функций введем искусственную вязкость q, определяемую следующим выражением q = PoI£kk(ciI£kk +c2C),если £kk > 0, q = 0, если £kk < 0,

где

I = V 13 , C =

Г 4G +д > ) ²

1 3 p o ep)

скорость звука в среде, c ₁ , c ₂ – безразмерные

константы, причем c 1 « 1,5, c ₂ =0,5-0,8.

В качестве определяющих соотношений, описывающих поведение рассматриваемых материалов, используем теорию течения [3].

Потеря сплошности среды происходит при достижении эквивалентным напряжением максимального значения.

Для решения краевой задачи (5)–(9), сформулированной в локальной форме, воспользуемся вариационным методом Галеркина [4]. Тогда искомое решение краевой задачи определяется из условия:

u ₀ = inf П ( u ), V t e D t , для { u e W^D ), u |_g D = ф } ,

где, П ( u ) есть слабая форма уравнения движения системы (5) с краевыми условиями (6), (7), (7а) и (8) на границе 5 D , u - вектор перемещений движущейся среды (деформируемого тела) c неоднородными граничными условиями.

б

г

в

Рис. 1 а) строение глаза, б) расчетная схема, в) конечно-элементная модель осколка, г) конечно-элементная модель глазного яблока: 1 – роговая оболочка, 2 – передняя камера, 3 – хрусталик, 4 – радужка, 5 – конъюктива, 6 – цилиарное тело с мышечными отростками и волокнами ресничного пояса, поддерживающее хрусталик, 7 – склера, 8 – сосудистая оболочка, 10 – сетчатка, 11 – стекловидное тело, 12 – центральная ямка желтого пятна, 13 – артерии, 14 – зрительный нерв, 15 – оболочки зрительного нерва, 16 – внутренняя прямая мышца глаза

Таблица 1

Вид повреждения	Скорость осколка, м/с	Характерный размер осколка L, мм
Разрыв роговицы	37	1,12
	58	1,8
Разрыв склеры	32	0,98
	63	0,69

σ , МПa

Рис. 2. Схематизированные диаграммы деформирования: 1 - роговицы, 2 - склеры

Рис. 3а. Отскок осколка от роговой оболочки

Задача (10) эквивалентна задачи отыскания минимума функционала

δΠ= ∫ ( pv i ^t -σ ij , j V

- pf_i ) δ u_idV + ∫ ( σ _ijn_j S_F

- t_i ) δ u_ids +

+ ∫ ( σ ¹ ij S_C

- σ ² _ij ) n _j δ u_ids = 0.

Для этого расчетная область D покрывается конечно-элементной сеткой D^h , составленной из конечных элементов Δ k таким образом, что:

M

D « Dh = J Ak, k=1

а в качестве пробных функций используются полиномы Лагранжа, определяемые на каждом конечном элементе в виде:

uj(xi) = Nr(δjk)uˆkr , где uˆkr – узловые значения искомого поля перемещений uj(xi) , Nr(δjk) – так называемые функции формы [5], δjk – символы Кронекера.

Роговая оболочка

Рис. 3б. Проникновение осколка в глазное яблоко

Роговая оболочка

Внутриглазная жидкость

Подставляя соотношение (12) в функционал (11) и минимизируя невязку, получаем полудискретное уравнение метод конечных элементов c нелинейной, зависящей от времени и решения правой частью:

[ M ] V = F B - ^ J B T 0 ( 8 , 8 ) dV e , (13) e

Ve где [M ] - матрица масс в общем случае размерностью 3N х 3N (N - число узлов конечноэлементной сетки), которая после распределения масс конечных элементов по узлам и при использовании линейных функций формы в (12) является диагональной, FB – вектор внешних и массовых сил, B – матрица производных от функции формы, Ve – объем конечного элемента.

Для решения динамической контактной задачи (5) с граничными условиями (8) используется численный алгоритм “контактного взаимодействия”, разработанный ранее [6] и обобщенный на задачи ударной динамики.

Учет физической нелинейности деформируемой среды осуществляется на основе известного и достаточно хорошо апробированного итерационного метода “последовательного догружения с возвратом по радиусу” [7].

Для интегрирования по времени уравнения (13) используем центральную конечно-разностную схему [8]. Вычислительные особенности приведенной выше задачи приведены в [9].

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

На рис. 1 а-г представлены расчетная схема и конечно-элементные модели осколка и глаза. Диаграммы деформирования для роговицы и склеры [11] аппроксимируем в виде билинейных кривых (рис. 2).

Рис. 3а изображает картину отскока осколка от роговицы, не повреждая ее. На рис. 3б показано, как проникает осколок в глазное яблоко (фрагмент из общей картины). Повреждения склеры аналогичны. Фрагменты картин даны на вертикальной плоскости симметрии глазного яблока.

В таблице 1 приведены условия повреждения глазного яблока в зависимости от скорости и характерного размера осколка.

Выводы

Представленные здесь результаты исследований проводились в рамках программы РАН «человек-машина-среда». Они дают полезную информацию о возможных механизмах и условиях повреждения глаза и, таким образом, для обоснования разработки и надлежащего выбора специальных средств защиты от глазных травм при определенных технологических процессах (шлифование, операции, связанные с обработкой стекла, и т.п.), боевых действиях (взрывы снарядов и т.п.) и аварийных ситуациях. Полученные в результате вычислительного эксперимента предельные скорости и размеры осколка согласуются с данными натурных экспериментов (клиническими данными) [12, 13].

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