Диффузионная модель разрушения элементов опорно-двигательного аппарата человека Текст научной статьи по специальности «Математика»

УДК 531/534: [57+61]

ДИФФУЗИОННАЯ МОДЕЛЬ РАЗРУШЕНИЯ ЭЛЕМЕНТОВ ОПОРНО-ДВИГАТЕЛЬНОГО АППАРАТА ЧЕЛОВЕКА

А.В. Чигарев1, А.В. Борисов2

1 Белорусский национальный технический университет, Республика Беларусь, 220027, Минск, Проспект Независимости, 65, e-mail: chigarev@rambler.ru

2 Филиал Государственного образовательного учреждения высшего профессионального образования «Московский энергетический институт (Технический университет) в г. Смоленске», 214013, Смоленск, Энергетический пр., 1, e-mail: borisowandrej@yandex.ru

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

Ключевые слова: уравнение Фоккера-Планка-Колмогорова, функция плотности вероятности, микротрещина, разрушение, вероятность, время достижения границы, математическое ожидание.

Введение

С развитием математических моделей опорно-двигательного аппарата человека возникла проблема моделирования его прочности, надежности и долговечности. Как показывает практика, первостепенное внимание необходимо уделять зонам контакта, их усталостному износу и разрушению. Опорно-двигательный аппарат человека функционирует в условиях со случайным образом изменяющимися обстоятельствами внешней среды и, соответственно, нагрузками. Поэтому для моделирования наиболее адекватными являются статистические методы. В научной литературе имеются работы по определению усталости на металлических и неметаллических веществах неживой природы. Проблемы прочности в живой природе исследованы менее полно. Этому способствует практическая невозможность натурных экспериментов на живом объекте. Исходя из этого, необходимо создавать модели износа и разрушения биологических объектов на основе синтеза известных моделей для неживой природы и учета в них явлений, присущих живой материи [9]. Один из таких подходов заключается в рассмотрении статистической модели разрушения опорнодвигательного аппарата человека.

© Чигарев А.В., Борисов А.В., 2012

Чигарев Анатолий Власович, д.ф.-м.н., профессор, заведующий кафедрой теоретической механики, Минск

Борисов Андрей Валерьевич, к.т.н., доцент кафедры высшей математики, Смоленск

09806267

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

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

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

Постановка задачи определения функции плотности вероятности появления

И РАСПРОСТРАНЕНИЯ ТРЕЩИНЫ В ОДНОМЕРНОМ СЛУЧАЕ

Рассмотрим кинематику распространения трещины. Проведем анализ одномерного, двумерного и трехмерного случаев.

Вначале проанализируем одномерный случай. Здесь возможны аналитическое (для некоторых случаев) и численное решения.

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

Рассмотрим для определенности тазобедренный сустав человека (рис. 1). Трещина начинает свое движение из положения О. Предположим, что в одномерном случае движение трещины происходит вдоль направления ОМ.

На рисунке х0 - положение начала микротрещины в начальный момент времени t0, х0 = х(^); х - положение конца микротрещины к моменту времени t , х = х^) .

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

Во-вторых, на сустав действует нерегулярная, стохастическая составляющая нагрузки, возникающая в результате случайных перегрузок, например во время падений, занятий спортом и т.п. В суставе имеют место неоднородности структуры материала, другие микротрещины и т.п. Под действием второго фактора в одномерном случае изменяется только скорость роста микротрещины, а в дву- и трехмерном

Рис. 1. Постановка задачи распространения трещины в тазобедренном суставе человека в одномерном случае (кортикальный слой кости значительно увеличен)

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

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

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

В одномерном случае возможно следующее описание распространения трещины с использованием уравнения Фоккера-Планка-Колмогорова [4, 5]:

t t)]+1 ±_ [К 2 t ^ )], (1)

где ^(х, t) - одномерная функция плотности распределения значения координаты конца трещины к моменту времени ^ К1 (х, t) - коэффициент сноса в биологическом материале; К2 (х, t) - коэффициент диффузии в биологическом материале; t - время распространения трещины; х - координата конца трещины. Начальные и граничные условия задаются в виде 5-функций, сосредоточенных в заданных точках.

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

Аналитическое решение уравнения Фоккера-Планка-Колмогорова

ДЛЯ ОДНОМЕРНОГО СЛУЧАЯ

Рассмотрим аналитическое решение уравнения (1). Колмогоровым [4] был предложен следующий метод решения с введением новых независимых переменных.

t1 = ), х1 = у(х, t). (2)

Плотность вероятности с использованием подстановки (2) принимает вид

^ = ™\ф(\у(хх|. (3)

После подстановки новых переменных и преобразования [4] уравнение (1) принимает вид

17 >] ) = -~х К (х!, t1 V: (х!, ^ )]+ | [К 2 (х!, t1 V: (х!, t1 )] . (4)

д^ дх1 2 дх1

Подбирая функции ф^) и у(х, t), можно определить коэффициенты К1 (х1, t1) и К2 (х1, t1) так, чтобы из уравнения (4) можно было найти w1 (х1, t1). В работе Колмогорова [5] предложены такие функции:

К1 (х, t) = а(^ )х + Ъ(;), К 2 (х, t) = с(^), (5)

где а(^), Ъ(^), с($) - функции времени.

Выбираем функции ф^) и у(х, t) в виде

t t

-J a(t ' )dt ' t _J a(t ' )dt ' t ~ 2 J a(t ' )dt '

y(x,t)= xe t0 -Jb(t ” )e'° dt”, ф(>)= Jc(t ")e '0 dt". (б)

t0 t0

При таком выборе функций ф(>) и y(x, t), как показано в работе [5],

коэффициенты в новых координатах равны

K1 ta t1 )= 0, K 2 (xl, t1 )= 1. (7)

Тогда уравнение Фоккера-Планка-Колмогорова принимает вид

Xwl (xl, tl) = 1a 2 wl (xl, tl) (8)

Xt1 2 йх.2

Это уравнение параболического типа. Его решение известно из теории уравнений математической физики [1]

(x1 - xl0 ^

w'(x" )= ь ( іe ^). (9)

VMt - >1,)

По формуле (3) получаем окончательное решение.

Применим данный метод к решению задачи о распространении трещины в твердых тканях (суставе) человека.

Рассмотрим в нулевом приближении модель, в которой функции а (t), Ъ(^), с^) являются константами. Пусть а(^) = а, Ъ(^) = Ъ , с^) = с . Выбор линейных видов

зависимостей для коэффициентов сноса и диффузии объясняется следующими обстоятельствами.

Коэффициент сноса описывает детерминированную часть процесса. Это регулярное конвективное ламинарное движение. Под действием регулярного усилия трещина начинает расти прямолинейно, перпендикулярно усилию.

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

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

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

Вычисляя функции ф^) и у(х, t) по формулам (6), имеем

Используя формулы (2), (3), (9), получаем искомую функцию плотности вероятности координаты положения конца трещины:

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

Теперь на основании выражения (11) проведем численный анализ полученного решения в системе компьютерной математики «Мяйетайса» (лицензионная версия данной программы приобретена за счет средств гранта президента Российской Федерации № МК-2524.2008.1). Выберем произвольным образом значения констант: а = 0,7; Ъ = 0,8; с = 0,5. Тогда поверхность решения w(x, t) имеет следующий вид

Построим сечения данного процесса при фиксированных значениях аргументов (рис. 3).

При анализе графиков сечения процесса и исходной зависимости (рис. 2, 3, а) видно, что наиболее вероятно появление микротрещин в первой половине времени процесса.

По рис. 2, 3, б видно, что вероятность появления трещины зависит от ее длины. С ростом длины трещины вероятность ее появления уменьшается. Рост длины трещины имеет сопротивление, и, чтобы трещина росла, необходимо прикладывать дополнительное внешнее усилие.

у(х, t) = хе а{‘-к) - (1 - е а{‘-к ])Ъ/а, ф^) = (1 - е-2а{{-к ])с!2а .

(10)

(11)

(рис. 2).

Рис. 2. Зависимость функции плотности вероятности w от х и от t

а б

Рис. 3. Зависимость функции плотности вероятности при фиксированных значениях:

а - х = 1; б - t = 1

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

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

По рис. 2 и 3 видно, что в результате получаются зависимости, близкие по форме к нормальному закону распределения, который является наиболее распространенным в природе.

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

Численное решение уравнения Фоккера-Планка-Колмогорова

ДЛЯ ОДНОМЕРНОГО СЛУЧАЯ. ПРОБЛЕМА ПРАВИЛЬНОГО ОРТОДОНТИЧЕСКОГО

ПЕРЕМЕЩЕНИЯ

Чтобы проверить достоверность полученных результатов, проведено численное решение исходного уравнения Фоккера-Планка-Колмогорова (1) в системе компьютерной математики «Mathematical [3]. Интегрирование проведем на отрезке от 0 до 6 по обеим переменным. Начальные и граничные условия задаем в виде

w(x,0) = e~x , w(0,t)= e'* , w(6,t)= e~6 .

42k 42k 42k

В результате получаем следующий график (рис. 4).

Сравнивая рис. 2 и рис. 4, видим достаточно убедительное сходство аналитического решения и численного решения, что доказывает адекватность используемых методов. Различия при x ^ 6 связаны со спецификой численного решения - заданием граничного условия.

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

Установлено, что наиболее вероятны появление трещин небольшого размера в начальный период и рост размера трещин с увеличением времени процесса.

Перейдем к анализу реальной картины, наблюдаемой в суставе. Как выше указывалось, вид коэффициентов сноса и диффузии при исследовании живого организма установить невозможно. Так как проведенный анализ показал устойчивость решений к изменению их вида, ограничимся случаем, когда данные коэффициенты являются константами. Предположим, что коэффициенты сноса и диффузии являются средним значением для химических элементов, входящих в состав кости. В справочной литературе указывается, что приводимые данные по коэффициенту диффузии «справедливы в лучшем случае по порядку величины» [2]. Поэтому примем усредненное значение для коэффициента диффузии K2 = 10 9 м2/с. Коэффициент сноса

оценим как K1 = ^K2 /1 м/с, где t - время процесса.

Рис. 4. Результат численного решения уравнения Фоккера-Планка-Колмогорова

а б

Рис. 6. Сечения функции плотности вероятности при фиксированных значениях:

а - х = 0,706; б -1 = 0,5

Время жизни человека в системе СИ порядка 109 с, размер в суставе порядка 10-3 м. До этих значений необходимо вести численное интегрирование уравнения Фоккера-Планка-Колмогорова, что сопряжено со значительным накоплением ошибки при их решении. Поэтому нормируем указанные две величины, чтобы интегрирование вести в пределах от 0 до 1. После соответствующих пересчетов коэффициенты сноса и диффузии будут равны: К1 = 1,6-10~7, К2 = 6,7-10 5 м/с.

Тогда в результате расчетов для одномерного случая получаются следующие зависимости (рис. 5).

Также построим сечения поверхности на рис. 5 при фиксированном одном из аргументов (рис. 6).

По данным графикам видно, что с течением времени функция плотности вероятности возрастает в процессе жизни. Зависимость функции плотности вероятности от длины трещины такая: вначале она быстро возрастает, затем начинает убывать. Это значит, что в начале жизни человека идет инкубационный период накопления большого количества микротрещин малого размера, затем их количество стабилизируется, и идет медленным темпом накопление локальных микротрещин, и только к концу жизни маловероятно возникновение магистральной трещины. Наиболее вероятно обнаружить трещину длиной 0,05-0,10 от длины исходного сустава, которая соответствует скрытой фазе накопления микроповреждений.

Определение функции плотности вероятности появления и распространения

ТРЕЩИНЫ В ДВУМЕРНОМ И ТРЕХМЕРНОМ СЛУЧАЯХ

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

Многомерное уравнение Фоккера-Планка-Колмогорова, записанное применительно к плотности вероятности состояний Н (X, t), имеет вид [4]

= -Е йТ К (x, * Мх *)]+ 1 Е Е йЛ- К (x, * Мх *)], (12)

й* ,=1 йх, 2 ,= ;=! йх,йх;

где: х - случайный вектор одномерных случайных функций х1 (*), х2(*), ..., хм(*), х = х(*); К1 (х, *) - коэффициенты сноса; Кц (х, *) - коэффициенты диффузии.

Применительно к двухмерному случаю в используемой прямоугольной декартовой системе координат имеем: х1 (* ) = х(*), х2 (* ) = У(*), а в трехмерном случае х3 (*) = г(*). N = 2,3 соответственно для двумерного и трехмерного случаев, х, (*) -координаты положения конца трещины вдоль соответствующей оси координат к моменту времени *. К (х, *) = К (х, у, *), К1 (х, у, г, *), К (х, *) = К (х, у, *), К (х, у, г, *) -коэффициенты сноса и диффузии соответственно для двумерного и трехмерного случаев. Для удобства данные коэффициенты будем обозначать индексами с координатами следующим образом: индексу с номером , = 1 будет соответствовать индекс х, с номером , = 2 - индекс у и с номером , = 3 - индекс г.

В двумерном случае уравнение Фоккера-Планка-Колмогорова в введенных обозначениях имеет вид

йн( х, у,*)

1

+-

2

Йх[Кх (x, у,*) Н(x, y, * )]+йу [Ку (x, y, *) Н(x, у,1:)]

+

й2 й

[Кхх (x, у^) Н(x, y, * )>^-[Ку (x, y, *) Н(x, у,*)]+ (13)

йх йхйу

й2 й2

+йуйх [Кух(x, у,г ) н( х У, * )]+й)7 [Куу(x, У, *) н( x, у,г )]

Предположим однородность и изотропность тканей, образующих сустав. Тогда

КХ = Кх = Ку , К 2 = ^ ^ = Ку = ^ .

Используя описанную ранее методику численного решения дифференциальных уравнений в частных производных, получаем функцию плотности вероятности в виде интерполяционных функций. Решая уравнение Фоккера-Планка-Колмогорова (13) при тех же числовых значениях, аналогично рассмотренному выше одномерному случаю, получаем сечения процесса при одном из фиксированных аргументов, изображенные на рис. 7.

Так как в двумерном случае функция плотности вероятности является функцией трех аргументов н(х, у, *), для ее построения требуется четырехмерное пространство. Поэтому построим ее сечения при фиксированном одном из аргументов. Численные значения берутся такие же, что и в одномерном случае, все коэффициенты сноса и диффузии предполагаются равными между собой соответственно. Фиксированный аргумент в каждом случае равен 0,5.

По рисунку видно сходство сечений при фиксированных х и у. Это объясняется указанными выше предположениями о равенствах коэффициентов вдоль осей. Таким образом, в суставе человека разрушение равновероятно может идти в любом направлении. Как и в одномерном случае, наиболее вероятно возникновение микротрещин в диапазоне 0,05-0,10 от длины исходного сустава. С течением времени (рис. 7, в) функция плотности вероятности ведет себя без значительных изменений, кроме начальных и конечных значений.

а

б

Рис. 7. Сечения функции плотности вероятностей при фиксированном одном из аргументов

Таким образом, решена двумерная задача нахождения функции плотности вероятности распределения координаты конца микротрещины к моменту времени

Теперь рассмотрим уравнение Фоккера-Планка-Колмогорова, которое для трехмерного случая имеет вид

х, у, ^,Ґ) 3 г / \ \ —I 3 г / \ 1

—*— зХ К' х у ^{^ х y, *,ґ '-_зу ^Ку 'x, у ^{^ х y, ^ К ' х y, 2,ґ' ^ х y, 2,ґ '-!+

+2 {!? К (х’у ^ ^х’уг’ ^ЗУ2 К (х’ У’ ^ ^ У’ ^я+З? К (х’ У’ ^ ™(х У’^ )]+ (14)

+3^К (х’ У’ 2’? ) ^ х’ У’2’? )]+зУзх К (х’ У’ 2’? ) ^ х’ У’2’? )]+зх& К (х’ У’ 2’* ) ^х’у’ 2’* )]+

з2 з2 з2 ]

+з23х К (х’ У’ 2’^(х’ У’2’^ ]+ЗУ32 К (х’ У’2’^(х’ У’2’^ ]+з23У К (х’ У’2’^(х’ У’]}.

Решаем уравнение Фоккера-Планка-Колмогорова для трехмерного случая при тех же числовых значениях и аналогично рассмотренным выше одномерному и двумерному случаям в предположениях К1 = Кх = Ку = К2 ’

К2 = Кхх = КУУ = К22 = КхУ = КУх = Кх2 = Кх = К^ = К2У ’ получаем функЦиЮ плотности

вероятности в виде функции четырех аргументов ^(х’ У’ 2’ ’ для ее построения требуется пятимерное пространство. Поэтому построим ^(х’ У’ 2’ ^) при фиксированных двух аргументах (рис. 8). Все значения берутся такие же’ как и прежде. Фиксированные аргументы в каждом случае равны 0’5.

Остальные графики не приводятся’ поскольку они являются идентичными двум построенным: графику а - если фиксируются х и У’ а также х и г; б - если фиксируются

I и х или I и у.

0,3'

0,3

м

а

б

Рис. 8. Графики функции плотности вероятностей при фиксированных двух переменных

По графикам видно’ что функция плотности вероятности в целом ведет себя аналогично рассмотренным выше более простым случаям.

В данном примере также имеется сходство зависимостей’ объясняемое предположением об однородности и изотропности материала. Таким образом’ проанализирован трехмерный случай.

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

Определение вероятности достижения трещиной границ сустава

Задача нахождения вероятности достижения трещиной границ сустава является важной в практическом плане. Пока есть только много мелких повреждений’ локальных микротрещин в элементах опорно-двигательного аппарата человека’ он будет работать. Когда мелкие трещины сливаются в магистральную трещину’ становится возможным явление перколяции (просачивания жидкости). В этом случае нормальное функционирование уже становится невозможным. Рассмотрим задачу достижения магистральной трещиной границы на примере коленного сустава.

Найдем вначале вероятность достижения односторонней границы. Подобная задача важна в случае отсутствия какой-либо части ткани’ например’ вследствие операции по установке искусственного сустава. Важным является оценивание вероятности разрушения кости в процессе эксплуатации имплантата.

Вероятность Р(х0^) того’ что траектория процесса достигнет за время I границы а или Ь’ равна [4]

где х0 - начальное значение координаты в начальный момент времени ї = ї0; Q(x0,ї) -вероятность того, что трещина не достигнет границы за время ї.

Функция Q(x0,ї) определяется с помощью функции q(x0, х, ї) следующим образом

Функция q(x0’ х’ {) не является плотностью вероятности в обычном смысле’ так как она не удовлетворяет условию нормировки’ но удовлетворяет граничным условиям

(15)

(16)

а

q(x0, а, ї) = q(x0, Ь, ї) = 0 .

(17)

Рассмотрим случай достижения односторонней границы на расстоянии L Ф 0.

В этом случае формула (15) с учетом (16) принимает вид

Ь

Р(х0, ^ = 1 #(х0, х, t)dx . (18)

—ад

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

Фоккера-Планка-Колмогорова при выбранных должным образом начальных условиях.

Определим функцию ^(х0, х, ^ в виде линейной комбинации плотностей вероятности следующим образом:

#(х0, х, ^ = ^(х0, х, {) — Ы(хг, х, {), (19)

где х1 - дополнительный полюс, расположенный в точке зеркального отражения

основного полюса х0 от границы £, т.е. х1 = 2£ - х0; k выбирается исходя из условия

t

|[1г(х0, £, ^ - kw(x1, £, t)]dt = 0. (20)

0

Тогда искомая функция д(х0, х, ^ после преобразований принимает вид

Ы\ е*1/2 -1^-ах + ах0еа‘2/2 ^ ь\е-“'^2-1^у/а-2L+х0 + хе

аІ2 ( аґ2 , Л аІ2 ( -аґ2 і

----- ае\ е -1 І - С е -1

/ \ е 2 е - е 2 ke

^(хо,х, І )=------------------------------------------------------------------------------------------------------. , л~--. (21)

^(1 - е аІ )/а

Проведем численный анализ полученного результата.

Находя численно значение k для полученной функции плотности вероятности при вышеуказанных значениях параметров и при размерах, соответствующих характерным размерам в суставе человека х0 = 0,01 м; £ = 0,02 м; х1 = 0,03 м, получаем k = 1,03.

Построим график функции д(х0, х, t) (рис. 9).

Вычисляя вероятность достижения трещиной границы (18) в течение жизни, получаем значение Р = 0,99.

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

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

Наличие фундаментального решения уравнения Фоккера-Планка-Колмогорова позволяет воспользоваться методом отражений [4].

Граничные условия в этом случае имеют вид

#(х0, £/ 2, t )= #(х0,-£/ 2, t )= 0. (22)

Рис. 9. Вид функции д ( хо, х, t) в зависимости от х

-0,2 -0,1 0,1 0,2 Рис. 10. Вид функции д(Хо,х,t) в зависимости от х

Функция д(хо, х, і) определяется следующим образом:

п=ад

д(хо, х, і )= ^(- 1)п^(п£, х, і).

(23)

Выписывать в явном виде полученную нами функцию д(хо, х, і) не будем ввиду

её громоздкости. Приведем её график на рис. 10.

В этом случае вероятность достижения трещиной границы принимает вид

Ц 2

Р(хо, і) = 1 - Гд(хо, х,і)dx.

-V 2

(24)

Вычисляя численно интеграл при тех же параметрах, что и ранее, получаем следующее значение вероятности: Р = 0,295. Это означает, что в процессе жизни человека вероятность разрушения элементов опорно-двигательного аппарата не очень велика.

Рассмотренный выше метод возможен в случае, если уравнение Фоккера-Планка-Колмогорова имеет фундаментальное аналитическое решение, доступное для анализа, т.е. только для одномерного случая. Во всех остальных случаях необходимо прибегать к численным методам решения.

п=-ад

Рассмотрим численное решение данной задачи. Так как численное решение уравнения Фоккера-Планка-Колмогорова получено, определим вероятность достижения трещиной размеров, равных размерам материала в суставе за время жизни человека. Как известно [6], такая вероятность равна

хт

Р(х, t) = Цw{x, t)dxdt . (25)

0 0

Вычисляя численно двойной интеграл в пересчитанных в относительных единицах значениях, получаем следующее значение: Р(х, t) = 0,297. Как видно, отличие от аналитически найденной вероятности минимально, менее 1%.

Оценим теперь численно вероятности достижения трещиной размеров сустава в двумерном и трехмерном случаях.

В двумерном случае

ХУТ

р(х, у, t ) = л/ ^(х, t^хёу&. (26)

000

Получаем значение: Р(х, у, t )= 0,360.

В трехмерном случае

ХУ Z т

Р(х, у, z, t ) = //// ^(х, t)dxdуdzdt. (27)

0000

Получаем значение: Р(х, у, z, t ) = 0,385.

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

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

Определение среднего времени первого достижения трещиной границ сустава

Найдем среднее время первого достижения трещиной границ сустава. Уравнение для времени достижения границы процессом Т (а, х0, Ь) получено

в предположении о стационарности процесса и приводится в работах [4, 8]. Воспользуемся приводимыми в них результатами.

і _ V ( \ ЛТ 2Т

_ 1 _ ^ ТК2 (х0 ) , 2 .

ах0 2 ах0

Аналитическое решение данного уравнения приводится там же:

(28)

и у И

Т =

еф(г а

е-ф(у )ау | е ~ф(у ]йу +| |

К 2 (Г)

е ф( г ]йг

и

е ~ф( у )ау | е~ф( у )ау

и

| е-ф( у )ау

х

0 1_

0 1_

Рис. 11. Зависимость среднего времени первого достижения трещиной границы от начального положения начальной микротрещины

Рис. 12. Зависимость среднего времени первого достижения трещиной границы от начального положения начальной микротрещины, численное решение

Проводя численные расчеты при значениях границ, соответствующих характерным размерам в суставе а = 0, Ь = 6, для времени достижения границы из соответствующего начального положения имеем следующую зависимость (рис. 11).

Проведем численную проверку полученного результата. Решим численно исходное дифференциальное уравнение (28) при тех же значениях параметров и нормальных граничных условиях. В качестве результата строим график (рис. 12).

По графикам видно хорошее совпадение результатов аналитического и численного решений.

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

Выводы

Построено аналитическое решение уравнения Фоккера-Планка-Колмогорова применительно к одномерной задаче распространения трещины в биологическом материале.

Проведено численное решение уравнения Фоккера-Планка-Колмогорова применительно к той же задаче, показана идентичность решений и возможность использования численных методов, реализованных в среде системы компьютерной математики, к решению задач распространения микротрещин.

Численно решены задачи распространения микротрещин в двумерном и трехмерном случаях.

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

Полученная функция плотности вероятности применена к оценке вероятности достижения трещиной границ тканей сустава, несущих нагрузку, в течение всей жизни человека. Установлено, что вероятность находится в пределах от 0,297 до 0,385.

Оценено математическое ожидание размеров микротрещины и времени ее появления. Установлено, что с переходом от одномерного случая к двумерному и трехмерному их значения возрастают, в целом находясь в пределах от 0,13 до 0,19 от исходного размера сустава и времени их образования в пределах от 0,15 до 0,20 от времени жизни человека.

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

1. Араманович И.Г., Левин В.И. Уравнения математической физики. - М.: Наука, 1969. - 288 с.

2. Григорьева И.С., Мейлихова Е.З. Физические величины: справочник. - М.: Энергоатомиздат, 1991. -1232 с.

3. Дьяконов В.П. Mathematica 5/6/7: полное руководство. - М.: ДМК-Пресс, 2008. - 624 с.

4. Казаков В.А. Введение в теорию марковских процессов и некоторые радиотехнические задачи. -М.: Советское радио, 1973. - 232 с.

5. Колмогоров А.Н. Аналитические методы теории вероятностей // Успехи математических наук. -1938. - Вып. V. - С. 283.

6. Кремер Н.Ш. Теория вероятностей и математическая статистика: учеб. для вузов. - М.: ЮНИТИ-ДАНА, 2007. - 573 с.

7. Миллер Б.М., Панков А.Р. Теория случайных процессов в примерах и задачах. - М.: Физматлит, 2002. - 320 с.

8. Тихонов В.И. Выбросы случайных процессов. - М.: Наука, 1970.

9. Чигарев А.В., Борисов А.В. Предельные нагрузки в суставах человека // Вестник ЧГПУ им. И.Я. Яковлева. Серия: Механика предельного состояния. - 2010. - № 2 (8). - С. 548-552.

A DIFFUSION MODEL OF ELEMENTS OF THE HUMAN LOCOMOTOR SYSTEM DESTRUCTION

A.V. Chigarev (Minsk, Belarus), A.V. Borisov (Smolensk, Russia)

A diffusion model of the destruction of elements of the human locomotor system is proposed in the work. Methods of the probability theory and mathematical statistics are applied to the description of the mechanisms of destruction of biological tissues. Propagation of microcracks is expected to be a Markovian process. Using diffusion equation of Fokker-Planck-Kolmogorov, the probability density function of distribution of coordinates of a microcrack end in hard tissues of the human musculoskeletal system is obtained. Practically, important statistical characteristics are defined on the basis of the above-mentioned function: the probability of achieving the joint boundary, the average time of the first reach of the boundary, mathematical expectation and the dispersion for coordinate values of a microcrack end and time of achieving them. One-, two- and three-dimensional cases are considered. Numerical estimations and their comparison with the real objects have been carried out.

Key words: Fokker-Planck-Kolmogorov equation, probability density function, microcrack, destruction, probability, time of reaching the boundary, mathematical expectation.

Получено 26 ноября 2011