CC BY
28
УДК 531/534: [57+61]
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ПОВРЕЖДЕНИЯ ГЛАЗНОГО ЯБЛОКА ПРИ КОНТАКТЕ С ЛЕТЯЩИМ ОСКОЛКОМ
А.И. Надарейшвили, В.А. Петушков, К.В. Фролов
Институт машиноведения РАН, Россия, 101990, Москва, М. Харитоньевский пер., 4, e-mail: aleko@mail.com
Аннотация. На основе уравнений ударной динамики, определяющих соотношений механики сплошных сред, с использованием конечно-элементных и конечноразностных аппроксимаций выполнен вычислительный эксперимент по проникновению осколка в глазное яблоко и определению механических и физических условий повреждения роговицы и склеры.
Ключевые слова: вычислительный эксперимент, модели, проникновение, метод конечных элементов.
Введение
Несмотря на неутешительные прогнозы о возможности восстановления нормального функционирования травмированного глаза, важными являются анализ причин и выявление механизмов повреждения глазного яблока с целью их дальнейшего предотвращения.
Причины, вызывающие повреждения глаза, могут иметь механическую, тепловую, химическую и т.п. природу. Большая доля глазных травм (до 70%) происходит на производстве [1]. Это указывает не только на уязвимость глаза, но и на недостаточные меры безопасности, применяемые против таких травм.
Среди производственных глазных травм наиболее опасным являются травмы с внутриглазным проникновением инородного тела. Такие повреждения глазного яблока обычно имеют летальный исход.
Многочисленные клинические данные о последствиях глазных травм, размеры извлекаемых осколков, а также анализ причин и условия возникновения летящих осколков можно считать результатами натурных экспериментов по проникновению инородного тела в глазное яблоко. На основе этих данных с использованием современных методов математического моделирования и компьютерной техники, появляется возможность изучения не только физических и механических условий глазных травм, но и оценки влияния этих травм на последующее функционирование поврежденного глаза.
В приведенном ниже материале моделируются процессы деформирования и разрыва тканей роговой оболочки и склеры глазного яблока человека при контакте с летящим инородным телом. Моделирование процесса проникновения осколка ограничивается разрывом роговицы (или склеры) и началом вхождением осколка в жидкостную среду первичной камеры.
© Надарейшвили А.И., Петушков В.А., Фролов К.В., 2005
В основе математической постановки задачи лежит подробный учет реального строения глазного яблока и свойств заполняющих его биологических веществ. Численное решение соответствующей краевой задачи строится с использованием комбинированного метода конечных элементов и конечных разностей соответственно для описания геометрии глаза и решения по пространству и на временном слое.
Постановка задачи. Математические модели
Летящее инородное тело (осколок) будем рассматривать в виде абсолютно твердой частицы клиновидной формы. Скорость полета осколка может определяться, например, как произведение угловой скорости вращения ротора инструмента на его радиус, или скоростью воздушного потока, переносящего абразивные частицы.
Глазное яблоко рассматривается в виде сферических оболочек для роговицы и склеры по границе, заполненных жидкообразной средой с внутриглазным давлением 20 мм рт. ст, снаружи она нагружена атмосферным давлением, действующим на роговицу и склеру.
Отметим, что в такой объемной постановке удается учесть все основные составные части глазного яблока и сложную трехмерную картину распространения в нем ударных волн и процессов повреждения.
Для пространственной аппроксимации передней камеры, хрусталика и стекловидного тела были применены пространственные восьмиузловые, а для роговицы и склеры четырехузловые безмоментные конечные элементы.
Представленное таким образом глазное яблоко рассматривается далее как нелинейно деформируемая среда объема V, занимающая в Я3 произвольную, в общем случае многосвязную область Б, ограниченную поверхностью £. В этом случае область Б включает в себя поверхности раздела в виде соединения разнородных материалов (сред) и возможных дефектов в виде отслоений и повреждений.
Начальную конфигурацию к0 области Б отнесем к прямоугольной декартовой системе координат хг, г = 1,2,3 (эйлеровы координаты). Движение (деформацию) тела Б относительно начальной конфигурации к0(Б) в любой произвольный момент времени г > 0 можно описать следующим отображением:
X- = С( хк, г),
X = (X1,X2,X3) е к0(Б), г > 0, к0 :Б ^Я3,
где Хк - лагранжевы координаты, а С - непрерывно обращаемое отображение из к0 в Я3.
Мерой деформации тела Б служит градиент деформации Е = дх/дХ, который согласно известной теореме о полярном разложении в теории конечных деформаций может быть представлен в виде
Е = Яи,
где Я - положительно определенный ортогональный тензор вращения, и -симметричный тензор растяжения.
Градиент скорости частицы тела Б равен:
Ь = ^ = д^дХ = ЕТЕ-', (1)
дх дХ дх
этот тензор может быть представлен как сумма тензоров скоростей вращения Ж и растяжения У:
Ь = Ж + У,
где
1 т
У = -(Ь + Ьт), (2)
Скорость тензора деформации Грина-Сен Венана Е = 0,5(и2 -1) определяется из соотношения:
Е = ЕтУЕ. (4)
Среда движется (деформируется) под действием импульсных ударных воздействий, в общем случае, задаваемых или в виде поверхностных сил Ец (х1, г),
х1 е ^ £Е , г е Б( = (0, т), т - время действия импульса, или начального распределения
скоростей у. (х1, 0) на части поверхности , причем £ = £Е и и £Е п = 0 . Здесь
и далее в отношении индексов и операций над индексами используются соглашения, принятые в тензорном исчислении.
Система уравнений импульсной динамики деформируемой среды, описывающая ее механическое поведение, может быть представлена в следующем виде [2]:
р ёй Е = Р0,
ё1у а + р / = ру1, на Б х Бг (5)
ег = 4 - (Р + У)¥г,
где р - плотность рассматриваемой среды, верхний знак 1 означает производную по
времени, /ц - массовые силы, у = и\ - вектор скорости, а, в - тензоры напряжений и
деформаций, соответственно, =агц - рдц - девиатор тензора напряжений, р = аи/3
- давление, е - полная энергия на единицу объема среды, V - относительный объем текущего состояния, q - объемная вязкость.
В начальный момент времени г = 0
Р = Р0, р = Р0, е = е„, х,(Хк,0) = Хк, х](Хк,0) = V,(Хк,0), (6)
где V, - начальное распределение скоростей.
Приведем граничные условия для деформируемой среды Б. Запишем их для известной формы начальной границы £0 в момент времени 1 = 10, £0 = £(10) :
х,(Хк, 1) = г,(1), Хк е £0, (7)
а1}п} = Е (Хк, 1), Хк е £0, (7а)
При наличии внутренней границы контактного разрыва е £0 к приведенным условиям (7) добавляются условия контактной непрерывности:
(аЦ -а'2-)пц = 0, Хк = Х2, Хк е £с0 хБ,, (8)
где цифрами 1, 2 обозначены тела (или разнородные материалы), находящиеся по обе стороны от границы £с.
Если при этом возникает возможность относительного движения указанных тел с границами £1 и £(°2, первоначально совпадающими с £°, условию отрыва общих точек поверхности в зоне контакта соответствуют нулевые усилия в них,
(Ц =а1 пц = 0).
Условию скольжения с трением соответствует (-дЕ^’' 1< Е? <уЕ1/ ), где
индексами 1 и п отмечены касательное и нормальное направления к поверхностям в
12 12 12 12 точке контакта, д - коэффициент трения, Еп’ =ац п,пц , Е(’ =ац п,1ц, 1ц - единичный
вектор касательной.
Для уравнения состояния деформируемой среды мы ограничимся так называемым квазиакустическим приближением, когда давление является только функцией объема, а ударная адиабата совпадает с изоэнтропой расширения [2]:
р=* (V - *)•
где К - объемный модуль упругости
Для сглаживания разрывов в профилях изменения плотности, напряжений и других полевых функций введем искусственную вязкость q, определяемую следующим выражением
q = РоIskk (ciIskk + c2CX если > 0
q = 0, если skk < 0,
(9)
где I = VV3, C =
скорость звука в среде, cV, c2 - безразмерные
V 3Ро 8Pj константы, причем c1 « 1,5, c2=0,5-0,8.
В качестве определяющих соотношений, описывающих поведение рассматриваемых материалов, используем теорию течения [3].
Потеря сплошности среды происходит при достижении эквивалентным напряжением максимального значения.
Для решения краевой задачи (5)-(9), сформулированной в локальной форме, воспользуемся вариационным методом Галеркина [4]. Тогда искомое решение краевой задачи определяется из условия:
u0 = inf П(u), Vt е Dt, для ju е W2v(D), u|8D = ф|, (10)
где, П(и) есть слабая форма уравнения движения системы (5) с краевыми условиями (6), (7), (7а) и (8) на границе 8D, и - вектор перемещений движущейся среды (деформируемого тела) с неоднородными граничными условиями.
Рис. 1 а) строение глаза, б) расчетная схема, в) конечно-элементная модель осколка, г) конечно-элементная модель глазного яблока: 1 - роговая оболочка, 2 - передняя камера, 3 -хрусталик, 4 - радужка, 5 - конъюктива, 6 - цилиарное тело с мышечными отростками и волокнами ресничного пояса, поддерживающее хрусталик, 7 - склера, 8 - сосудистая оболочка, 10 - сетчатка, 11 - стекловидное тело, 12 - центральная ямка желтого пятна, 13 -артерии, 14 - зрительный нерв, 15 - оболочки зрительного нерва, 16 - внутренняя прямая
мышца глаза
Таблица 1
в
г
Вид повреждения Скорость осколка, м/с Характерный размер осколка Ц, мм
Разрыв роговицы 37 1,12
58 1,8
Разрыв склеры 32 0,98
63 0,69
а,
МПа
Рис. 2. Схематизированные диаграммы деформирования: 1 - роговицы, 2 - склеры
Рис. 3а. Отскок осколка от роговой оболочки Задача (10) эквивалентна задачи отыскания минимума функционала
5П = |(ру\ - ац,ц - р£ )Ьи^ + | (ацпу -1 )8чЖ +
V £е
+1 (а1^ - а2ц)пц5и ^ = 0.
(11)
Для этого расчетная область Б покрывается конечно-элементной сеткой Бк, составленной из конечных элементов Ак таким образом, что:
м
Б * Бк = У Ак,
к=1
а в качестве пробных функций используются полиномы Лагранжа, определяемые на каждом конечном элементе в виде:
иц(X) = Nг (5 Цк) Ы к,
(12)
где и к - узловые значения искомого поля перемещений Ыц (х), Ыг (5 цк) - так называемые функции формы [5], 5цк - символы Кронекера.
£
С
Внутриглазная жидкость
Рис. 3б. Проникновение осколка в глазное яблоко
Подставляя соотношение (12) в функционал (11) и минимизируя невязку, получаем полудискретное уравнение метод конечных элементов с нелинейной, зависящей от времени и решения правой частью:
[М]V- = Ев-£{Вта(г,е)Л',, (13)
где [М] - матрица масс в общем случае размерностью 3Ы х 3Ы (И - число узлов
конечноэлементной сетки), которая после распределения масс конечных элементов по узлам и при использовании линейных функций формы в (12) является диагональной, Ев - вектор внешних и массовых сил, В - матрица производных от функции формы,
V - объем конечного элемента.
Для решения динамической контактной задачи (5) с граничными условиями (8) используется численный алгоритм “контактного взаимодействия”, разработанный ранее [6] и обобщенный на задачи ударной динамики.
Учет физической нелинейности деформируемой среды осуществляется на основе известного и достаточно хорошо апробированного итерационного метода “последовательного догружения с возвратом по радиусу” [7].
Для интегрирования по времени уравнения (13) используем центральную конечно-разностную схему [8]. Вычислительные особенности приведенной выше задачи приведены в [9].
Результаты моделирования
В результате вычислительного эксперимента установлены размеры и скорости проникновения осколка через роговую оболочку и склеру глазного яблоко. Плотности тканей роговицы и склеры равны 1149 и 1243 кг/м3 [10], а для внутриглазной жидкости передней камеры и стекловидного тела 1000 кг/м3. Хрусталик глаза рассматривается как абсолютно жесткое тело.
На рис. 1 а-г представлены расчетная схема и конечно-элементные модели осколка и глаза. Диаграммы деформирования для роговицы и склеры [11] аппроксимируем в виде билинейных кривых (рис. 2).
Рис. 3 а изображает картину отскока осколка от роговицы, не повреждая ее. На рис. 3б показано, как проникает осколок в глазное яблоко (фрагмент из общей картины). Повреждения склеры аналогичны. Фрагменты картин даны на вертикальной плоскости симметрии глазного яблока.
В таблице 1 приведены условия повреждения глазного яблока в зависимости от скорости и характерного размера осколка.
Выводы
Представленные здесь результаты исследований проводились в рамках программы РАН «человек-машина-среда». Они дают полезную информацию о возможных механизмах и условиях повреждения глаза и, таким образом, для обоснования разработки и надлежащего выбора специальных средств защиты от глазных травм при определенных технологических процессах (шлифование, операции, связанные с обработкой стекла, и т.п.), боевых действиях (взрывы снарядов и т.п.) и аварийных ситуациях. Полученные в результате вычислительного эксперимента предельные скорости и размеры осколка согласуются с данными натурных экспериментов (клиническими данными) [12, 13].
Надо отметить, что выполненные разработки для исследования процессов деформирования и разрушения в сплошных биологических средах могут быть успешно применены для изучения аналогичных процессов, протекающих в других тканях человеческого организма, таких как костные ткани и соединения при прыжках, падениях и пулевых воздействиях, в том числе при наличии бронезащиты и др.
Список литературы
1. Johnston, S.S. Symposium on ocular trauma / S.S. Johnston // Trans. Ophthalmol. Soc. - UK. - 1975. - 95.
2. Зельдович, Я.Б. Физика ударных волн и высокотемпературных гидродинамических явлений / Я.Б. Зельдович, Ю.П. Райзер. - М.: Физматгиз, 1966.
3. Малинин, Н.Н. Прикладная теория пластичности и ползучести / Н.Н. Малинин. - М.: Машиностроение, 1975.
4. Флетчер, К. Численные методы на основе метода Галеркина / К. Флетчер. - М: Мир, 1988.
5. Стренг, Г. Теория метода конечных элементов / Г. Стренг, Дж. Фикс. - М.: Мир, 1977.
6. Петушков, В.А. Моделирование нелинейных процессов разрушения повреждаемых упругопластических сред при высокоскоростном ударном нагружении двухфазной газожидкостной средой / В.А. Петушков, А.Н. Мельситов // Математическое моделирование. - 2000. - Т. 12. - № 10.
7. Helleux, J.P. Non-linear transient analysis for thin axisymmetric structures / J.P. Helleux // Numerical methods for non-linear problems. - 1980. - V. 1.
8. Рихтмайер, Р. Разностные методы решения краевых задач / Р. Рихтмайер, К. Мортон. - М., 1972.
9. Петушков, В.А. Математическое моделирование нелинейных процессов деформирования и разрушения объемных тел при высокоскоростном соударении / В.А. Петушков, А.И. Надарейшвили // Математическое моделирование. - 2004. - Т. 16. - N. 11.
10. Hogan, M.J. Histology of the human eye / M.J. Hogan, J.A. Alvarado, J.E. Weddell. - Philadelphia: Saunders, 1971.
11. Buzard, K.A. Introduction of biomechanics of cornea / K.A. Buzard //Refract Corneal Surg. - 1992. - 8.
12. Schrader, W. Perforating injuries: causes and risks are changing: a retrospective study / W. Schrader // Ger. J. Ophtalmol. - 1993. - V. 2. - P. 6-82.
13. Liggett, P.R. Ocular trauma in an urban population / P.R. Liggett, K.J. Pince, W. Barlow // Ophthalmology. - 1990. - V. 97. - P. 581-584.
MATHEMATICAL MODELING OF DAMAGE OF THE EYEBALL AT CONTACT WITH THE FLYING SPLINTER
A.I. Nadareishvili, V.A. Petushkov, K.V. Frolov (Moscow, Russia)
On a basis of the equations of the shock dynamics, constitutive relations of continuum mechanics, with use of finite element and finite difference approximations computational experiment with penetration of a splinter in the eyeball and with definition of mechanical and physical conditions of damage of the cornea and sclera is executed.
Key words: computational experiment, models, penetration, finite element method.
Получено 10 июня 2005
