Full Text

ВВЕДЕНИЕ В современных исследованиях [17-18] подтверждается наличие предварительных напряжений (ПН) [6, 19] во многих биологических тканях. Причинами возникновения ПН могут быть как ростовые деформации, происходящие, например, в костной или хрящевой тканях, так и последствия травм и хирургических операций в мягких тканях. В качестве примера естественных предварительных © Богачев И.В., Недин Р.Д., 2021 Богачев Иван Викторович, к.ф.-м.н., научный сотрудник Института математики, механики и компьютерных наук имени И.И. Воровича, Ростов-на-Дону Недин Ростислав Дмитриевич, к.ф.-м.н., доцент Института математики, механики и компьютерных наук напряжений в биологических тканях можно рассматривать линии Лангера, являющиеся линиями естественного натяжения кожи [14]. Их наличие объясняется структурой кожи, содержащей ненатянутые коллагеновые волокна, которые при натяжении деформируются и частично распрямляются, обеспечивая подвижность кожи, например, при движениях мышц и конечностей. Современные экспериментальные данные по исследованию деформации кожного покрова [14] позволяют сделать вывод, что его поведение соответствует вязкоупругому композиционному материалу, состоящему из неоднородных слоев эпидермиса, дермы и подкожного жира, характеристики которых различаются. В работах [2, 10] подробно описана методика практического измерения упругих модулей кожного покрова в результате проведения экспериментов по одноосному и двуосному растяжению образцов кожного покрова, а также по релаксации и ползучести. Собраны данные для различных участков кожного покрова, измерения проводились для образцов кожи пациентов разного возраста. В работе [1] подробно описан и экспериментально проиллюстрирован процесс релаксации кожного покрова при больших деформациях. Статья [15] посвящена, в том числе, исследованию влияния предварительных деформаций на восстановление кожного покрова после операций. Что касается ПН, возникающих в кожном покрове после операционных вмешательств, то в качестве примера можно привести операцию лифтинга, рассмотренную в работе [13]. В рамках лифтинга, то есть подтяжки кожи, происходит сепарирование лоскута кожного покрова, который затем натягивается в одном направлении, и избыток кожи иссекается, после чего разрез ушивается. При этом в процессе релаксации остаточные (предварительные) напряжения выходят на определенный уровень, который необходимо контролировать на различных стадиях заживления, чтобы не допустить омертвения кожи вследствие нарушения кровотока (при чрезмерном натяжении) либо провисания (при недостаточном). Таким образом, отметим, что для качественной диагностики состояния кожного покрова, в частности, в послеоперационный период необходима разработка уточненных моделей, учитывающих неоднородность, реологические свойства и возможное наличие ПН. При этом для исследования живой кожной ткани у пациента необходимо использование неинвазивных методов, например методов акустического (ультразвукового) зондирования. Целями данной работы являются: · исследование деформаций кожного покрова, моделируемого вязкоупругим неоднородным слоем, в свою очередь состоящим из трех слоев, с учетом ПН; · разработка методики идентификации механических характеристик кожи и предварительных напряжений на основе неинвазивного акустического подхода. ПОСТАНОВКА ЗАДАЧИ Для вывода постановки задачи о колебаниях предварительно напряженного вязкоупругого слоя, моделирующего кожный покров, в данной работе будем использовать линеаризованную постановку задачи об установившихся колебаниях предварительно напряженно-деформированного анизотропного тела [6]. Ввиду того, что тело нагружено периодически меняющейся нагрузкой qe-iwt , уравнения колебаний и определяющие соотношения после отделения временного множителя примут вид [4, 7]: T + rw2u = 0; (1) ij , j i T = s + u0 s + u s0 . (2) ij ij i,m mj i,m mj mj В формуле (2) функции s0 описывают предварительные напряжения; величины без индекса - компоненты соответствующих добавочных векторов и тензоров: ui - компоненты вектора перемещений, sij - компоненты симметричного тензора напряжений Кирхгофа; Tij - компоненты несимметричного добавочного тензора напряжений Пиолы; r - плотность. Будем рассматривать кожный покров как изотропный слой, занимающий область S = {x1 , x3 Î(-¥, ¥), x2 Î[0, H]} в условиях плоской деформации. Тогда у вектора смещений ненулевыми являются компоненты u1 =u1 (x1, x2 , w) , u2 =u2 (x1, x2 , w) . Исходя из реальной конфигурации кожного покрова будем считать, что исходный слой состоит из трех составных слоев, толщины которых известны, и их механические характеристики являются функциями поперечной координаты x2 . Также ввиду того, что наиболее характерным в рассматриваемом случае является одноосное сжимающие поле ПН, определяющие соотношения (2) можно записать в форме T = l(u + u ) + 2mu + s0u ; 11 1,1 2,2 1,1 1,1 T22 = l(u1,1 + u2,2 ) + 2mu2,2 ; (3) T = m(u + u ), T = m(u + u ) + s0u . Функции l(x2 ), m(x2 ) 12 1,2 2,1 21 1,2 2,1 2,1 характеризуют законы изменения параметров Ламе. В качестве граничных условий будем рассматривать модельные условия жесткого защемления нижней границы слоя, как прикрепления кожного покрова, и нагружения нормальной нагрузкой q = (q1(x1), q2 (x1), 0) на верхней границе: 2 2 u1 x =0 = u2 x =0 = 0; 2 T12 x =H = q1 (x1 ), T22 x2 =H = q2 (x1 ). (4) Для описания вязкоупругих свойств кожного покрова используем модель стандартного вязкоупругого тела [9]. Применяя принцип соответствия [16] к определяющим соотношениям теории упругости (3), функции параметров Ламе l(x2 ) и m(x ) заменим на функции комплексных модулей l* (x , iw) = inwl2 (x2 ) + l1 (x2 ) , 2 2 1+ inw m*(x , iw) = inwm2 (x2 ) + m1 (x2 ) , зависящих от частоты колебаний. Введенные здесь 2 1+ inw величины: l1(x2 ),m1(x2 ) и l2 (x2 ),m2 (x2 ) - длительные и мгновенные модули соответственно, n - время релаксации, являются параметрами, описывающими модель. Постановка обратной задачи [5] заключается в идентификации механических характеристик кожного покрова (функций длительных и мгновенных модулей l1(x2 ),m1(x2 ) и l2 (x2 ),m2 (x2 ) , входящих в состав комплексных аналогов параметров Ламе) и предварительного напряжения по дополнительной информации о полях смещений, измеренных на поверхности слоя в заданном частотном диапазоне: 2 u1 x =H = f1 (x1 , w), u2 x =H = f2 (x1 , w), wÎ[w1, w2 ]. (5) 2 При этом форма, вид и характер зондирующего нагружения, описываемого граничными условиями (4), в эксперименте может варьироваться. Ввиду специфики рассматриваемой задачи для кожного покрова плотность r и параметр времени релаксации n считались постоянными и известными величинами. ПРЕДВАРИТЕЛЬНЫЕ НАПРЯЖЕНИЯ В ВЯЗКОУПРУГОМ ТЕЛЕ Для моделирования предварительного напряжения рассмотрим нестационарную задачу об определении напряженно-деформированного состояния кожного покрова в начальном состоянии (до операции, без ПН). Будем рассматривать постановку, характерную для пластической операции лифтинга, когда реализуется одноосное растяжение сепарированного лоскута, после которого фиксируется деформация натянутого лоскута кожи (с удалением излишков и ушиванием) и затем происходит процесс релаксации [13]. В соответствии с моделью стандартного вязкоупругого тела [9] определяющие соотношения, связывающие функции напряжений и деформаций, имеют вид E n ¶e + E e = n ¶s + s. (6) 1 ¶t 2 ¶t При этом в рассматриваемом случае, как уже отмечалось, величины мгновенного и длительного модулей E1 = E1 (x2 ) и E2 = E2 (x2 ) являются функциями поперечной координаты. С функциями l1(x2 ),m1(x2 ) и l2 (x2 ),m2 (x2 ) они связываются по обычным формулам для изотропного тела: Ei = mi (3li + 2mi ) (li + mi ) . При моделировании процесса релаксации при фиксированной деформации e = e0 соотношение (6) примет вид n ¶s + s = E e . (7) ¶t 2 0 Начальным условием в данном случае будет значение растягивающей нагрузки s t =0 = s1. Общее решение задачи (7), (8) имеет вид (8) s(x , t) = E2 (x2 )e0 + e-nt æs - E2 (x2 )e0 ö . (9) 2 n ç 1 n ÷ è ø Затем в рассматриваемой ситуации по истечении достаточно длительного промежутка времени, t ®¥ , напряжение в соответствии с решением (9) релаксирует до уровня, определяемого длительным модулем величиной начальной деформации e0 : E2 (x2 ) , временем релаксации n и lims(x , t) = E2 (x2 )e0 = s0 (x ) , (10) t ®¥ 2 n 2 значение которого затем используется в соотношении (3) для слоя с ПН. ПОСТРОЕНИЕ СХЕМЫ РЕШЕНИЯ ЗАДАЧИ Для решения задачи воспользуемся преобразованием Фурье [11] по переменной x1 , общий вид которого представлен ниже: ¥ Z (a, x2 ) = ò Z (x , x ) e-iax1 dx . 1 2 1 -¥ (11) Запишем уравнения задачи (1)-(3) в каноническом виде в трансформантах: T = ail*T (l* + 2m*) + (a2 (4m*(l* + m*) (l* + 2m*) - s0 )- rw2 )u ; 12,2 22 1 (12) Граничные условия (4) в трансформантах примут вид (13) Проведем обезразмеривание задачи (12)-(13) с помощью введения следующих обозначений для параметров и переменных: x2 = Hxˆ2 , x1 = Hxˆ1 , a = Haˆ , u1 = iHU1 , u = HU , T = im Tˆ , T =m Tˆ , l* = m l , m* = m m , s0 = m sˆ 0 , r = r rˆ , m = m (H ) , 2 2 12 0 12 22 0 22 0 0 0 0 0 1 r = r(H ) , q (x ) = HP(xˆ ) , w2 = m k2 H 2r : 0 i 1 i 1 0 0 Tˆ = aˆ lˆTˆ (l + 2m) + (aˆ 2 (4m(l + m) (l + 2m) - sˆ 0 )- rˆk2 )U ; 12,2 22 1 Tˆ = -aˆ Tˆ - (rˆw2 - aˆ 2sˆ 0 )U ; 22,2 12 2 (14) Далее для удобства записи знак "^" будем опускать. Представим функции, входящие в состав задачи (14), согласно методике, предложенной в [6], в виде разложений по параметру преобразования Фурье a , общий вид которых записан ниже: Затем выпишем задачи при различных степенях a : (15) (16) (17) a2 : T (2) = l*T 12,2 22 1 1 (1) (l* + 2m*) + (4m*(l* + m*) (l* + 2m*) - s0 )U (0) - rk2U (2); T (2) = -T (1) + s0U (0) - rk2U (2) ; 22,2 12 2 2 U (2) = T (2) m* + U (1) ; (18) 1,2 12 2 U (2) = T (2) (l* + 2m*) -l*U (1) (l* + 2m*); 2,2 22 1 U (2) 1 x2 =0 = U (2) x2 =0 = 0, T12 (2) x2 =1 = P(2) , T (2) x2 =1 = P (2) . 2 1 22 2 Рассмотрим задачи (16) при a0 . Заметим, что в них не входит функция ПН s0 , однако относительно функций (0) U i и Ti 2 (0) они являются несвязанными, и для них можно построить схему последовательного определения функций l* и m* , решая обратные задачи, подобные исследованным ранее в работе [3]. Для решения удобно последовательно рассмотреть два случая точечного нагружения: P1(x1) = Pd(x1), P2 (x1) = 0; P1(x1) = 0, P2 (x1) = Pd(x1), (19) (20) где d(x1) - дельта-функция Дирака, при которой трансформанты функций нагружения не будут зависеть от параметра преобразования, что упрощает процедуру разложения по параметру a . В таком случае для решения обратных задач используется получаемая из условия (5) дополнительная информация вида Затем с применением метода (21) линеаризации формулируется итерационный процесс, состоящий из решения прямых задач методом пристрелки [8] для текущих приближений искомых функций длительных и мгновенных модулей и последующего решения систем интегральных уравнений Фредгольма 1-го рода относительно поправок к искомым функциям с помощью метода регуляризации Тихонова [12]. Итерационный процесс продолжается до выполнения условия малости величины невязки. Детальное описание схемы итерационного процесса приведено в [3]. Далее рассмотрим задачу (18) при a2 в случае распределенной нагрузки ìïconst, x1 Î[-c, c]; [ ] P1 (x1 ) = 0, P2 (x1 ) = í ïî0, x1 Ï -c, c . (22) Нагрузка такого рода (нормальная к поверхности) является наиболее применимой для практической реализации. С другой стороны, ввиду четности, ее трансформанта Фурье также является четной вещественной функцией, соответственно, в разложении (15) для такой нагрузки ненулевыми будут только коэффициенты при четных степенях a . Следовательно, величины с индексом "(1)", являющиеся множителями при a1 , обращаются в нуль. Тогда задача (18), так же как и задача (15), 2 2 разделяется на две несвязанные задачи, при этом одна из них имеет нулевые граничные условия, а другая, записанная относительно трансформант функций перемещения, имеет вид ((l* + 2m*)U 2,2 (2) ) ,2 = s0U (0) - rk2U (2) ; (23) Условие (5) о дополнительной информации для обратной задачи в данном случае примет форму: i s U (2) = 2 x2 =1 f (2) (k ), ks Î[k- , k+ ]. (24) Заметим, что в задаче (23) функцию (0) U 2 необходимо определять из решения прямой задачи (16) для нагрузки вида (22). Также обратная задача определения функции ПН s0 после подстановки соотношения (10) в уравнение (23) сводится к определению вещественного параметра начальной деформации e0 . Для этого воспользуемся модификацией разработанного ранее проекционного метода алгебраизации [20]. Выпишем слабую постановку задачи, введя функцию V2 - поле возможных перемещений, скалярно домножив на нее уравнения (23) и воспользовавшись граничными условиями 1 1 1 E e (l* + 2m*)U (2)V dx - k2 U (2)V dx + 2 0 U (0)V dx - P (2)V (1) = 0. (25) ò 2,2 2,2 2 ò 2 2 2 ò n 2 2 2 2 2 0 0 0 Введем в рассмотрение jn - систему линейно независимых функций, удовлетворяющих однородным граничным условиям задачи (23). Решение краевой задачи (23) будем искать в виде U (2) (x , k) = a + åa (k)j (x ), (26) N 2 2 0 i i 2 i=1 где комплексные коэффициенты ai ( i = 1..N ) зависят от частоты колебаний k , слагаемое a0 обеспечивает выполнение неоднородного граничного условия задачи (23) о нагружении на поверхности. При подстановке представления (26) в соотношение (25) и выборе в качестве функции V2 = jn ввиду линейности его левой части по переменным U2 и V2 получим систему N линейных уравнений относительно коэффициентов которой представлен ниже: ai (k) , общий вид N åa A (k) = P (2)j (1) - e 1 E 2 U (0)j dx , n = 1..N. (27) i in 2 n 0 ò n 2 n 2 i=0 0 Затем полученная линейная система (27) решается относительно коэффициентов ai (k) . Таким образом, из разложения (26) получается выражение для U (2) (x , k, e ) , 2 2 0 зависящее от частоты k и искомого параметра e0 , которое затем подставляется в дополнительное условие (24) обратной задачи. Выбирая некоторое значение km , из условия (24) получаем уравнение для определения e0 : 2 m 0 U (2) (1, k , e ) = f (k m ). (28) Подводя итог, заметим, что идентификация всех искомых величин может проводиться в два этапа: 1. Идентификация функции комплексного модуля m* и последующая идентификация суммы модулей l* + 2m* , по которой определяется комплексный модуль l* из задач (16), с использованием функций нагружений (19), (20) и дополнительной информации (21). Для решения используется итерационный процесс, позволяющий восстанавливать в том числе функции, имеющие разрывы, что характерно для рассматриваемой модели трехслойного кожного покрова. 2. Идентификация функции s0 с использованием найденных на первом этапе функций комплексных модулей l*,m* или суммы l* + 2m* . Задача сведена к определению вещественной величины начальной деформации, возникшей в результате эксперимента на релаксацию e0 . Идентификация основана на использовании зондирующей нагрузки вида (22) и дополнительной информации (24). Для решения адаптирован специальный проекционный метод. ВЫЧИСЛИТЕЛЬНЫЕ ЭКСПЕРИМЕНТЫ Как отмечалось в статьях [13, 14], средний уровень модулей кожного покрова составляет ~ 5×106 Па, а кровоток в коже прекращается при нормальном напряжении max s = 104 Па [13]. Операция по лифтингу должна проводиться таким образом, чтобы был достигнут максимальный косметический эффект, при этом напряжения релаксировали до величины, меньшей smax . В связи с этим в экспериментах обезразмеренный параметр уровня ПН в кожном покрове после операции можно считать величиной, находящейся в диапазоне s0 ~ 10-2,10-3 Па. Проведем анализ влияния уровня ПН, определяемого начальной деформацией e0 , на акустический отклик - функции амплитудно-частотных характеристик (АЧХ) на поверхности слоя. На рис. 1 представлены графики АЧХ для различных значений уровня ПН: e0 = 0 (отсутствие ПН, сплошная линия) и 0 e = 10-3, 5 ×10-3,10-2 (пунктирные). Из графиков видно, что различия в АЧХ имеются, однако существенными они становятся лишь вблизи максимумов АЧХ (соответствующих резонансам в упругом случае). В связи с этим в качестве частотного диапазона зондирования на втором этапе (как и на первом, что проиллюстрировано в [3]) следует выбирать частоты вблизи вязкоупругого резонанса. Рис. 1. Влияние уровня ПН на АЧХ для различных значений e0 Далее представлены результаты идентификации функций длительных и мгновенных модулей для комплексных аналогов параметров Ламе и последующей идентификации начальной деформации, определяющей предварительные напряжения. Для трехслойного кожного покрова функции длительных и мгновенных модулей считались кусочно-линейными функциями: ì0,175 + 0, 5x, x £ 0, 5895; m (x) = ï0, 5 + 0, 2x, x £ 0,8655; ì0,19 + 0, 5x, x £ 0, 5895; m (x) = ï0, 55 + 0, 2x, x £ 0,8655; 1 í 2 î ï0,8 + 0, 2x, x £ 1; ì0, 2 + 0,8x, x £ 0, 5895; í î ï0,85 + 0, 2x, x £ 1; ì0, 3 + 0,8x, x £ 0, 5895; l (x) = ï + 0, 5x, x £ 0,8655; l (x) = ï + 0, 5x, x £ 0,8655; 1 í1,1 î ï1,8 + 0, 4x, x £ 1; 2 í1, 2 î ï1, 9 + 0, 4x, x £ 1. На первом этапе идентификация самих функций мгновенных и длительных модулей проводилась в диапазоне между первым и вторым резонансами, как наиболее информативном для идентификации, что показали результаты других проведенных экспериментов. Для восстановления функций mi потребовалось 8 итераций, погрешность не превосходила 3 % (рис. 2), для функций li + 2mi - 11 итераций, погрешность также не превосходила 3 %. Итоговая погрешность для вычисленных функций li не превосходит 6 % (рис. 3). На рис. 2 и 3 графики искомых функций представлены сплошной линией, графики восстановленных - точками. На втором этапе на основе полученных на первом этапе значений комплексных модулей Ламе проводилась идентификация параметра уровня ПН для различных его значений. Точность идентификации сравнивалась для двух различных значений частот, на которых в эксперименте снималась дополнительная информация (24): вдали от первого резонанса ( k= 0.1) и возле него ( k =1.23). Результаты, представленные в таблице, показывают, что точность идентификации существенно выше для частоты, близкой к вязкоупругому резонансу. В скобках указана погрешность идентификации. Рис. 2. Этап 1. Результаты восстановления функций mi Рис. 3. Этап 1. Результаты восстановления функций l i Идентификация параметра уровня ПН № Точное e0 k= 0,1, (погр., %) k = 1, 23, (погр., %) 1 0,05 0,04961 (0,77) 0,04983 (0,35) 2 0,01 0,00974 (2,62) 0,00988 (1,15) 3 0,005 0,00475 (4,94) 0,00489 (2,17) 4 0,001 0,00077 (23,4) 0,00089 (10,3) Также отметим, что при уменьшении величины уровня ПН снижается его влияние на АЧХ, в связи с чем ухудшается точность его определения. Точность также зависит от погрешности идентификации комплексных модулей на первом этапе. ЗАКЛЮЧЕНИЕ Разработана двухэтапная методика определения механических характеристик кожного покрова и параметра уровня одноосного предварительного напряжения. Методика позволяет определять комплексные модули с погрешностью, не превосходящей ~ 6 %, что является удовлетворительным результатом для подобного рода задач. Также она показала свою применимость для восстановления уровня ПН в диапазоне, соответствующем значениям релаксировавших ПН в кожном покрове после операций, на практике. Наибольшая точность восстановления происходит при зондировании на частотах, близких к резонансным.

About the authors

I. V Bogachev

R. D Nedin

References

Statistics

Views

Abstract - 133

PDF (Russian) - 65

Refbacks

• There are currently no refbacks.