CC BY
21
УДК 539.3
О КОЛЕБАНИЯХ ФУНКЦИОНАЛЬНО-ГРАДИЕНТНОЙ ПОРОУПРУГОЙ КОЛОННЫ*
© 2014 г. А.О. Ватульян, А.А. Ляпин, Ю.А. Святко
Ватульян Александр Ованесович - доктор физико-математических наук, профессор, заведующий кафедрой теории упругости, факультет математики, механики и компьютерных наук, Южный федеральный университет, ул. Мильчакова, 8а, г. Ростов н/Д, 344090; заведующий отделом дифференциальных уравнений, Южный математический институт Владикавказского научного центра РАН, ул. Маркуса, 22, г. Владикавказ, 362027, e-mail: vatulyan@math.sfedu.ru.
Vatulyan Aleksandr Ovanesovich - Doctor of Physical and Mathematical Science, Professor, Head of Department of the Elasticity Theory, Faculty of Mathematics, Mechanics and Computer Sciences, Southern Federal University, Milchakov St., 8a, Rostov-on-Don, 344090, Russia; Head of Differential Equation Division, Southern Institute of Mathematics of Vladikavkaz Scientific Center RAS, Marcus St., 22, Vladikavkaz, 362027, Russia, e-mail:
vatulyan@math.sfedu.ru.
Ляпин Александр Александрович - кандидат физико-математических наук, младший научный сотрудник, Научно-исследовательский институт механики и прикладной математики Южного федерального университета, пр. Стачки, 200/1, г. Ростов н/Д, 344090, e-mail: lyapin@sfedu.ru.
Lyapin Alexander Alexandrovich - Candidate of Physical and Mathematical Science, Junior Scientific Researcher, Research Institute of Mechanics and Applied Mathematics of Southern Federal University, Stachki Ave, 200/1, Rostov-on-Don, 344090, Russia, e-mail: lyapin@sfedu.ru.
Святко Юлия Александровна - магистр, кафедра теории упругости, факультет математики, механики и компьютерных наук, Южный федеральный университет, ул. Мильчакова, 8а, г. Ростов н/Д, 344090, e-mail: yuliannasvyatko@mail. ru.
Svyatko Yulia Alexandrovna - Master Student, Department of the Elasticity Theory, Faculty of Mathematics, Mechanics and Computer Sciences, Southern Federal University, Milchakov St., 8a, Rostov-on-Don, 344090, Russia, e-mail: yuliannasvyatko@mail. ru.
Изучена задача о колебаниях функционально-градиентной пороупругой колонны в режиме установившихся колебаний. Показана достаточная точность построения решения задачи методом пристрелки и проекционным методом типа Галеркина. При помощи линеаризации исходной системы дифференциальных уравнений сформулированы интегральные соотношения для решения обратной задачи по восстановлению переменного коэффициента проницаемости среды. Проведен ряд численных экспериментов по восстановлению.
Ключевые слова: пороупругость, обратные задачи, численные методы.
Using numerical methods for solving .systems of differential equations with variable coefficients the problem for steady-state vibrations of poroelastic functionally-gradient column is studied. For two numerical methods: shooting method and Galerkin kind method the corresponding accuracy is shown. Using linearization approach for base system of differential equations the corresponding integral equations for solving inverse problem for reconstruction ofpermeability coefficient are formulated. The group of numerical experiments for reconstruction of unknown parameters are realized.
Keywords: poroelasticity, inverse problems, numerical methods.
Исследованию динамики сред сложной структуры в настоящее время уделяется все больше внимания. Одним из распространенных объектов исследования является пороупругая насыщенная жидкостью среда, позволяющая эффективно моделировать механиче-
ское поведение грунтов и биологических объектов, таких как костная ткань.
Одной из наиболее популярных моделей при описании динамики пороупругой среды является модель Био. В рамках данной модели уже выполнено доста-
*Работа выполнена при поддержке гранта РФФИ №14-01-31453 мол_а, а также в рамках программы фундаментальных исследований по стратегическим направлениям развития науки Президиума РАН № 1 «Фундаментальные проблемы математического моделирования».
точное количество работ, посвященных колебаниям одномерных тел (стержней, колонн) [1], легких конструкций из пористых материалов [2-5], исследованию процессов деформаций костных тканей, а также некоторых мягких тканей - хрящей, связок, сухожилий [6-10].
При моделировании сложных функционально-градиентных сред значительную трудность представляет задача определения материальных параметров. Зачастую при разработке очередной математической модели в механике деформируемого твердого тела, учитывающей более сложное поведение или новые свойства среды, возникают дополнительные параметры или функции, которые необходимо определять из соответствующих экспериментов. При отсутствии таковых одним из эффективных методов определения таких параметров является формулировка и решение обратной задачи, заключающейся в решении соответствующих краевых задач с учетом некоторой дополнительной информации.
Важной характеристикой пороупругой насыщенной жидкостью среды является параметр проницаемости. Задача об определении таких характеристик в зависимости от координат важна и актуальна.
Постановка задачи
Для описания движения данной среды в режиме установившихся колебаний рассмотрим модель Био, описывающую установившиеся колебания с частотой ю вертикальной пороупругой неоднородной колонны высотой к под действием приложенной к верхней грани нормальной нагрузки д. Данная модель является одномерным вариантом общей модели Био [11].
Краевая задача относительно амплитуд вертикальных смещений и и порового давления р в таком случае примет вид
[((Л + 2/и)ы ')'—ар'+ рю2и = 0,
\-irnap — (кр )'- 1ааи' = 0.
(1)
(2)
тему. Такое различие для некоторых пористых материалов может достигать порядка 1020^1030, что при решении соответствующих систем линейных уравнений будет приводить к значительным вычислительным погрешностям при решении жестких систем. При этом необходимо осуществить корректное обезразмеривание задачи.
Осуществим обезразмеривание системы уравнений (1) с граничными условиями (2). Введем следующие параметры:
и = к ■ V, х = к ■£, р = Р ■ Р*,
Л + 2 и = (Л + 2иШ) , ,
к (х) = ц(4)ко , ко = к (к),
2 рю2 к2
Лд + 8(£) =
aP
Лд + 2Мд
Л + 2м0 ф2 h
Л0 + 2 Ид рк2
v(£) =
\
рк2 R(x) к0 Лд+2И~д к2
a , J = ■
q
Рк к0Р Л + 2 и о
Здесь величины Л ,и - значения функций упругих модулей в некоторой произвольной точке колонны, обычно характеризующие максимальное значение функций Л,и; Р* - некоторое характерное давление в среде.
В результате обезразмеренную систему дифференциальных уравнений относительно V и Р можно представить в виде
(gV')' -yP' + k2V = 0, -(^P')' - iKÖP - IkvV' = 0. Граничные условия в безразмерном виде
V = 0, P = 0 при £ = 0, gV'-yP = J, P' = 0 при £= 1.
(3)
(4)
Граничные условия для задачи (1): [ и = 0, р = 0 при х = 0,
\(Л + 2и)и'—ар = д,р' = 0 при х = к.
Здесь Л, и, а, к - соответственно упругие модули, модуль Био и коэффициент проницаемости, которые в рамках исследуемой задачи являются функциями координаты х; д характеризует интенсивность нагрузки на верхнем торце колонны х=к.
В дальнейшем для проверки результатов будем решать поставленную задачу (1), (2) численно при помощи двух методов: пристрелки и Галеркина.
Основная проблема при численной реализации задач пороупругости - значительное различие в порядках материальных параметров, входящих в сис-
Описание методов решения краевой задачи
Первый метод, который будет реализован для решения краевой задачи (3), (4), - метод пристрелки. Для линейных систем дифференциальных уравнений его реализация заключается в решении промежуточных вспомогательных задач Коши для системы (3) и последующего представления решения в виде линейной комбинации решений промежуточных задач с некоторыми коэффициентами, причем эти коэффициенты линейной комбинации находятся из удовлетворения условию на верхнем торце колонны £ = 1.
Сведем задачу к канонической системе дифференциальных уравнений первого порядка, вводя новые неизвестные функции, характеризующие обобщенное напряжение и поток поровой жидкости в колонне.
(5)
S = gV-yP , Q = -VP', V ' = S-У P,
S ' = -k2V ,
P ' = -1Q, v
S У
Q = IKSP + IKV{ - + y P);
V(0) = 0 , P(0) = 0 , Q(l) = 0, S(1) = J . (6) Общее решение краевой задачи (5), (6) можно представить в виде линейной комбинации линейно независимых решений задач Коши, описываемых системой уравнений (5) и граничными условиями:
V1 (0) = 0, P(0) = 0, Qi(0) = 0, Si(0) = 1, (7)
V2(0) = 0, P2(0) = 0, Q2(0) = 1, S2(0) = 0. (8)
Задачи (5), (7) и (5), (8) решались численно с помощью метода Рунге - Кутты четвертого порядка.
В результате функции продольного смещения V и порового давления P можно представить в виде линейной комбинации V = AV + Лу2, P = A1P1 + A2P2. Коэффициенты Al, Л будем искать из условия на верхнем торце Q(1) = 0, S(1) = J , что приводит к простой системе линейных уравнений, которая решена методом Крамера.
Одним из наиболее эффективных методов решения краевых задач является проекционный метод типа Галеркина, который позволяет относительно просто строить решение в виде некоторой линейной комбинации базисных функций. Такое представление решения эффективно еще тем, что позволяет в соответствующих случаях получать не только решение задачи, но и строить производные от него, что в случае метода пристрелки можно реализовать только после соответствующей сплайн-аппроксимации решения.
Введем в рассмотрение дифференциальные операторы
L : (gV') '-yP' + k2V = 0 , Z2: -(vP') ' - iKÖP - IkvV' = 0.
(9)
(10)
Будем искать решения, удовлетворяющие граничным условиям
V(0) = 0, р(0) = 0,
8©V'(1) -уР(1) = J , Р'(1) = 0.
Построим решение задачи (9), (10) методом Галер-кина. В соответствии с теорией метода введем два набора базисных функций: фк, \ук. Тогда представления для функций смещения и давления можно записать в виде
N N
V = Х афф, +Ф0 , Р = £+^0 . (11)
,=1 ,=1
В граничных условиях (10) присутствует смешанное условие, характеризующее значение обобщенного напряжения на верхнем конце колонны. В таком слу-
чае автоматически удовлетворить этому условию не представляется возможным. Одним из способов преодоления трудности в этой ситуации является использование слабой постановки задачи. В качестве произвольных функций, участвующих в слабой постановке, будем использовать те же самые базисные функЦии фк , ¥к •
Умножим оба уравнения из системы (9) на функции ф, \ук соответственно. С учетом оставшихся граничных условий в (9) функции ц/к, фк выберем в
_ г Jd%
виде Г0
¥к = sin(A t).
0 g (t)
¥0 =0; Фк = sin<X 4),
Далее считаем, что /лк = \ =-^ + жк, к = 1,2...
С учетом граничных условий и после соответствующих преобразований получим систему
1
^' - уР)фк (1) -1 (8V 'ф'к - уР' ф[ + к^фк ^ =0, (12)
0
1
-^Р'¥к (1) + 1 (лРУк - 1к8Р¥к -'¥к)d# = 0 • (13)
0
Учитывая разложения функций V, Р по базисным ф , \ук, получим систему линейных алгебраических уравнений относительно коэффициентов разложения
а, и ъ ,
ХЛ а + Bl ь s = fk
k=1
£л1а, + Blbs = fk
(14)
где
Aks =-\W + ^ФА )dt + (gф'фk )| ,
0
1
Al =J (-iKV¥sWk)dt ,
0 1
B1 =J (-У¥фФ1 )d t- (У¥ А )| ,
0
1
Bl = J(V¥'s¥'k -ixÖYWk)dt-(V¥'s¥k)|t=1
fl = (Фo'Фk)
=1 -J (gФo' Ф'к +^ФфФк )dt,
fk = J(-1кф¥к)dt.
Решая систему (14), находим коэффициенты и Ь„ после чего - решение исходной задачи.
Осуществим ряд вычислительных экспериментов в соответствии с описанными выше методами. В табл. 1, 2 представлены результаты расчетов в зависимости от числа координатных функций N.
s=1
0
t
0
0
Таблица 1
Результаты для метода Галеркина в зависимости от числа N
Частота
Амплитуды к=0,1 к=2
N
3 5 9 3 5 9
0,787888 0,787886 0,787886 2,171338 2,156350 2,153341
lP®l 0,027367 0,027297 0,027284 0,830131 0,813365 0,809619
Нетрудно видеть, что при низких частотах для достижения точности 10-3 достаточно использовать всего лишь три координатные функции. Для более
высоких частот для получения приемлемой точности обычно достаточно 7-9 координатных функций, причем их количество растет при приближении к резонансным ситуациям.
На рис. 1 изображены графики функций смещения и давления при величине безразмерной частоты к=2, построенные различными методами для экспоненциальных зависимостей.
б
а
Рис. 1. Графики смещения (а) и порового давления (б):----- действительная часть, метод пристрелки;----- мнимая
часть, метод пристрелки;...........- действительная часть, метод Галеркина;---мнимая часть, метод Галеркина
На основе предложенных вычислительных схем построены и исследованы амплитудно-частотные характеристики (АЧХ) колонны. Выявлено наличие ограниченных резонансных пиков, амплитуды которых снижаются с ростом частоты. Данные об АЧХ далее использованы для решения обратных задач.
Таблица 2
Разница значений, полученных двумя методами
Здесь АV(1) = \^(1) — Vг(1)|, АР(1) = |р(1)—Р2(1)|; Р V соответствуют решению методом пристрелки, Р V - методом типа Галеркина.
Решение обратной задачи
Рассмотрим задачу о восстановлении коэффициента проницаемости среды ц как функции координат. Решение прямой задачи для неоднородной среды является неотъемлемым этапом при решении обратной задачи по восстановлению соответствующих коэффициентов дифференциального оператора [12]. Обратная задача является нелинейной и для ее решения используется итеративная процедура, основанная на методе линеаризации. Для этого необходимо составить интегральные уравнения относительно малых добавок к искомым коэффициентам, определяемых на каждом шаге итерационного процесса.
Будем считать, что в обратной задаче известна информация вида
V(1,к) = (к), Р(1,к) = /2(к), к^[кх,к2 ] (15) и что в краевой задаче (5), (6) функции g, у известны, а ^ - неизвестна.
Амплитуды Законы изменения
g(£) = 1+#2, v£) = 1+£ g (£) = e£ ,
к=0,1 к=1 к=0,1 к=1
т.1) 0,000728 0,012823 0,000407 0,006334
AP(X) 0,000125 0,009925 0,000077 0,005553
Проводим линеаризацию, делая замену щ = щ+ ещ, V = V , Р = Р0 +еР1, и формулируем
краевые задачи при степенях е0,е'.
е°: (V' -ГР' + ^К = 0, -ЩоР)' - К8Р0 - ^ = 0, К и=о, Ро1х=о=0, (8К-гРо)1=1 = 1, -Щ Ро'1 х=1 = 0, Пи=/1 К), Ро1х=1=/2 К). е1: ' -ГРХ1 = о,
-(щР) - ЩР') -гк5Рх -КУ? = о, пи=о, Р1 х=о=о, ( ¿V -гР)1==1=о, (-ЩР'-ЩР') и = о,
VI и=о, Р1\х=1=о.
Производя необходимые преобразования, используя граничные и дополнительное условия, исключая поправки V и Р , получим интегральное уравнение Фредгольма первого рода [13] для поправки:
Щ(4)[Р^] ^ = -/К, (16)
Ке[К,К2 ] ,
где /0 (к) - априорная информация о поровом давлении на верхнем торце колонны для функции щ , измеренная в соответствующем диапазоне частот [к, к2 ].
31
1
п 1,5
О
Заключение
Изучено динамическое поведение неоднородной пороупругой колонны. На основе данных об АЧХ, а также с помощью метода итерационной регуляризации построено решение обратной задачи о реконст-
Отметим, что ядро интегрального уравнения (16) обращается в нуль при 4 = 1, что приводит к плохой реконструкции искомой функции в окрестности этого края.
Общий подход к формулировке таких уравнений для задач пороупругости и восстановление других переменных характеристик описаны в работе [14].
Таким образом, решение обратной задачи сводится к последовательности решений интегральных уравнений Фредгольма 1 -го рода методом Тихонова с автоматическим выбором параметра регуляризации [15] в рамках итерационного процесса, который останавливается при выполнении определенного условия точности, наложенного на выражение разности АЧХ для эталонной задачи и текущей.
Приведем результаты ряда вычислительных экспериментов по восстановлению коэффициента проницаемости . В качестве входного значения данного параметра были использованы функции щ(4) = 4 + 4 + 1 (рис. 2а) и = 2£3 + 4£2 + 8£ + 3 (рис. 2б). Сплошной линией на графиках обозначено точное решение, пунктиром - начальное приближение (линейная функция). Важным параметром при решении обратной задачи является диапазон частот, на котором происходит съём априорной информации. Для представленных результатов были выбраны диапазоны [к1,к2\ = [8,71; 13,53] для рис. 2а и
[К, к2 ] = [4,91;16,91] для рис. 2б, расположенные между первым и вторым резонансами.
1
рукции переменного коэффициента проницаемости. Проведен ряд численных экспериментов по реконструкции неоднородного коэффициента проницаемости среды, показана достаточная точность при восстановлении.
■ у •' -у // // у/ •
.'• // //
-1 1-1-1 ---1 1 ' 1
0,2 0,4 0,6 0.8 1
4
б
Рис. 2. Результаты восстановления функции щ(х) :------ восстановленная характеристика;
--точная характеристика;..................- начальное приближение
а
Литература
1. Bardet J.P. The damping of saturated poroelastic soils during steady-state vibrations // Appl. Math. And Comp. 1995. Vol. 67. P. 3 - 31.
2. Allard J.-F. Propagation of sound in porous media: Modeling sound absorbing materials. London, 1993. 284 p.
3. Leclaire P., Horoshenkov K.V., Cummings A. Transverse vibrations of a thin rectangular porous plate saturated by a fluid // J. Sound and Vibrat. 2001. Vol. 247, № 1. P. 19 - 31.
4. Leclair P., Horoshenkov H.V., Swift M.J. The vibrational response of a clamped rectangular porous plate // J. Sound and Vibrat. 2001. Vol. 247, № 1. P. 19 - 31.
5. Theodorakopoulos D.D., Beskos D.E. Flexural vibrations of poroelastic plates // Acta Mech.1994. Vol. 103. P. 191 - 203.
6. Cowin S.C. Bone Poroelasticity // J. Biomech. 1999. Vol. 32, № 3. P. 217 - 238.
7. Suh J.K., DiSilvestro M.R. Biphasic poroviscoelastic behavior of hydrated biological soft tissue // J. Appl. Mech. 1999. Vol. 66. P. 149 - 156.
8. Маслов Л.Б., Блескин Е.В., Вихрев С.В. Исследование вынужденных колебаний биомеханической системы «кость-фиксатор» // Биомеханика-2000 : тез. докл. V Всерос. конф.
по биомеханике. 29 мая - 2 июня 2000 г., г. Н. Новгород. Н. Новгород, С. 34 - 35.
9. Маслов Л.Б., Ликсонов Д.В. Биомеханические характеристики нижней конечности человека // Математические модели и компьютерное моделирование в биомеханике : учеб. пособие / под ред. А.В. Зинковского и
B.А. Пальмова. СПб., 2004. С. 285 - 438.
10. Маслов Л.Б. Резонансные свойства большеберцовой кости в неповрежденном состоянии и с устройством внешней фиксации // Рос. журн. биомеханики. 2003. Т. 7, № 2. С. 20 - 34.
11. Biot M.A. General theory of three-dimensional consolidation // J. of Applied Physics. 1941. Vol. 12. P. 155 - 164.
12. Ватульян А.О. Обратные задачи в механике деформируемого твердого тела. М., 2007. 224 с.
13. Ватульян А.О. Интегральные уравнения в обратных задачах определения коэффициентов дифференциальных операторов теории упругости // Докл. РАН. 2005. Т. 405, № 3. С. 343 - 345.
14. Ватульян А.О., Ляпин А.А. Об обратных коэффициентных задачах пороупругости // Изв. РАН. МТТ. 2013. № 2.
C. 114 - 121.
15. Тихонов А.Н., Арсенин В.Я. Методы решения некорректных задач. М., 1986. 287 с.
Поступила в редакцию 9 июля 2014 г.
