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

Владикавказский математический журнал 2016, Том 18, Выпуск 2, С. 31-40

УДК 539.3, 539.5, 519.63

О ЗАДАЧЕ КОШИ В ТЕОРИИ КОЭФФИЦИЕНТНЫХ ОБРАТНЫХ ЗАДАЧ ДЛЯ УПРУГИХ ТЕЛ1

А. О. Ватульян, Л. С. Гукасян, Р. Д. Недин

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

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

Введение

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

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

© 2016 Ватульян А. О., Гукасян Л. С., Недин Р. Д.

1 Работа выполнена при поддержке Программы фундаментальных исследований по стратегическим направлениям развития науки Президиума РАН № 1 «Фундаментальные проблемы математического моделирования» (114072870112), «Математическое моделирование неоднородных и многофазных структур», Гранта Президента Российской Федерации МК-5440.2016.1, и Российского фонда фундаментальных исследований, проект № 16-38-60157 мол_а_дк.

границе при известном виде воздействия (вторая постановка) [1]. В обоих случаях необходимо решать некоторые обратные задачи на основании подходов, изложенных в работе [2, 3]. В первом случае обратная задача по определению неизвестных функций (обычно это коэффициенты Ламе) линейна, во втором — существенно нелинейна.

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

1. Общая постановка обратной задачи для упругой среды

Рассмотрим установившиеся колебания упругой среды в двумерном случае. Пусть плоская область Б С И,2 ограничена гладкой кривой I = 1\ и ¡2- Краевая задача об установившихся колебаниях с частотой ш имеет вид

(сцы пк,1)з + рш2щ = 0, г = 1,2, (1)

Пг1 к =0, Сцк1 пк,1 Щ 12 = Рг, (2)

где Сцк1 (х) — компоненты тензора упругих постоянных, удовлетворяющие известным свойствам симметрии, р(х) — плотность, которые являются функция ми координат XI, Х2, иг — компоненты вектора смещений1. Задача (1)-(2) о нахождении компонент вектора иг достаточно подробно исследована в литературе при известных гладких и кусочно-гладких законах неоднородности, и для нее доказаны теоремы существования и единственности [4], разработаны эффективные вычислительные схемы, основанные либо на конечноэлементных, либо конечноразностных подходах. В то же время в рамках этой модели весьма актуальными являются обратные задачи об определении законов изменения материальных характеристик по некоторой дополнительной информации о решении [5]. При этом наиболее часто исследуются две постановки. В первой постановке известны поля смещений внутри области на некоторой частоте, во второй — поля смещений на границе в некотором частотном диапазоне. В рамках первой и второй постановок исследовано достаточно большое число задач для упругих балок и стержней, когда колебания описываются обыкновенными дифференциальными уравнениями второго или четвертого порядка, однако для существенно двумерных и трехмерных задач такие задачи даже в первой постановке практически не исследованы [6].

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

2. Постановка обратной задачи в изотропном случае

Планарные колебания область Б описываются следующей краевой задачей:

огзз + рш2иг = 0, г, ] = 1,2, (3)

1 Индекс, стоящий после запятой, означает производную по соответствующей координате; например, ик,1 следует понимать как дик/дх\.

Ы\к =0, О;Щ \12 = Рг, (4)

где О; = \Ukk5ij + р(и»; + П;,г).

Здесь в постановку задачи входят две положительные функции А(Ж1,Ж2) и р(ж1 , Ж2), характеризующие законы изменения коэффициентов Ламе.

Основное внимание в настоящей работе уделим задаче определения модулей А(ж1,ж2) и р(жьж2) по известным компонентам поля перемещений в наборе точек из задачи 2. В развернутом виде краевая задача (3)-(4) имеет вид

(А(П1,1 + П2,2)),1 +(^(П1,2 + П2,1)) ,2 + (р(2иМ )),1 +р^2П1 = 0, (А(П1,1 + П2,2)),2 +(^(П1,2 + П 2,1)) ,1 + (р(2П2,2 )),2 +Р^П2 = 0,

П4 =0, (5)

[(А(П1,1 + П2,2) + р(2П1,1 ))П1 + р(П1,2 + П2,1 ^ = РЬ [(А(П1,1 + П2,2) + р(2П2,2))П2 + р(П1,2 + ^д)^] = Р2.

А (ж1, ж2)

р (Ж1, Ж2), удовлетворяющих (5) при известныX функциях П1 (Ж1,Ж2) И П2 (Ж1, Ж2) •

3. Исследование задачи о планарных колебаниях

Обезразмерим уравнения и граничные условия в (5), вводя безразмерные параметры и переменные

р А 2 Ро^2а2 р

5 =—, д = —, к, =-, ро = шах/л, а = шаш6, ро = тахр, 7 =—.

ро ро ро хе^ ро

Далее будем рассматривать случай, когда

9,5 е С~(£), С+°(5) = {/ : / £ С~(£), / ^ /о > 0}.

Тогда краевая задача (5) примет вид

0119.1 +Й12 5,1 +0135,2 +Ф1(ж,9,5) = 0,

0219.2 +Й225,1 +0235,2 +^2(ж, 5) = 0,

р; = р, р-1, (6) [(90ц + 5012)п1 + 5й1зп^ \12 = Р*, [(90ц + 5023)п2 + 5«1зп^ \12 = Р2*, где введены следующие обозначения

011 = П1,1 + П2,2 = 021, 012 = 2П1,1, 013 = П1,2 + П2,1, 022 = 013 , 023 = 2 П2,2,

Ф1 = 7К2П1 + 9011,1 + 5(012,1 + 013,2), (7)

Ф2 = 7К2П2 + 9021,2 + 5 (022,1 + 023,2).

Отметим, что система (6) представляет собой краевую задачу для системы дифференциальных уравнений в частных производных первого порядка относительно двух неизвестных функций q, g. Для того чтобы сформулировать задачу Коши для этих функций, необходимо, исходя из (6), определить значения искомых функций на границе ¿2- Для их определения необходимо, чтобы система (6) имела положительные решения; для разрешимости ее необходимо и достаточно, чтобы

А = an [ai3 (n2 - n|) + (ai3 - ai2) П1П2] = 0, (8)

а положительность регулируется ограничениями на нагрузку P*, P|. Так, например,

условие положительности нарушается, когда П2 = 0 P2 = 0 либо щ = 0, P* = 0, что

g

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

Если данные задачи таковы, что можно определить положительные данные Коши для искомых функций, то задача Коши

aiiq,i +ai2 g,i +ai3 g,2 +$i(x,q,g) = 0,

a2iq,2 +a22g,i +a23g,2 +$2(x, q, g) = 0, (9)

qln=o = qo, g|ra=o = go,

является определенной. В этом случае изучим вопрос о ее единственности. Для этого в (9) перейдем к локальной системе координат, связанной с частью границы области So С S, где вектор нагрузок (P* , P2) отличен от пуля. Пусть уравнение границы задано параметрически xi = xi(s), ж2 = x2(s), s £ [S_,S+]. Введем локальную систему координат (s,n) соотношениями

xi = xi(s) — nx'2(s), ж2 = x2(s) + nxi(s), (10)

считая, что (xi)2 + (ж^)2 = 1 и перепишем (9) в описанной системе координат. Тогда имеем следующую систему

Ацр + А12р + А13р- + Аир + Ф*(ж, q,g),

dn ds dn ds . .

(11)

A21P + A22p + A23p + A24p + Щх, q,g). dn ds dn ds

Разрешив эту систему относительно нормальных производных ^ и получим каноническую систему (при условии AiiA.23 — A2iAi3 = 0) следующего вида

dq , / dq dg \ dg , / dq dg \

^ = 4*'^ J ' ^ = ^ Г""to' ~ds) ' (12)

q|n=0 = qo' g|n=0 = go •

В предположении бесконечной дифференцируемости функций q, g решение зада-

n

ственность решения [7, 8].

Этот результат справедлив для некоторой подобласти So С S, которая получается, если использовать замену переменных (10). Для оставшихся подобластей S можно ис-

So

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

4. Некоторые аспекты численной реализации

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

Наиболее эффективными методами решения обратных задач в рамках сформулированной постановки являются метод конечных элементов и различные вариации проекционных методов [6]. Ниже, в п. 4.1, рассмотрим проекционный метод, основанный на использовании слабой постановки исходной краевой задачи (3), (4). Для этого введем гладкие функции V, удовлетворяющие главным граничным условиям, налагаемым на вектор перемещений: ^^ = 0. Умножая уравнения (3) па пробные функции V и интегрируя по области Б, получим

J {?V + рш2щV) ^ = 0.

я

Преобразуя интеграл с помощью формулы Грина, получим соотношение

J (о] V + рш2щ^) ^ = J О] ЩП] ^ - J (о] ^ - рш2им)

Я дЯ Я

которое с учетом граничного условия (4) примет вид

J РМ ^2 - J (о] V] - рш2и^) = 0. (13)

¡2 Я

Отметим, что способ получения слабой постановки исходной краевой задачи показывает, что для функции щ, удовлетворяющей главным граничным условиям, можно требовать гораздо меньших ограничений на гладкость, чем для функций сильного решения (т. е. решения задачи в ее исходной постановке) [9, 10], что важно при использовании конечноэлементных и проекционных трактовок.

4.1. Исследование прямой задачи. Прямая задача заключается в нахождении функций перемещения щ(ж1;ж2), и2(ж1,ж2) с учетом заданных коэффициентов Ламе. Воспользуемся идеей проекционных методов и представим искомые функции в виде конечного разложения по некоторым базисным функциям, удовлетворяющим главным граничным условиям. В качестве базиса возьмем набор бигармонических полиномов. Таким образом,

п

щ = ^ (14)

¿=1

где ф1 = Ж1Ж2, ф2 = ж1ж2, фз = Ж1Ж2, = Ж1 ж|, Ф5 = Ж4Ж2, фб = Ж^ ж|, ф/ = Ж^,

ф8 = Ж1 Ж5, фд = ЖбЖ2, ф10 = ж1 ж2, . . . ,

п

и2 =^2 , (15)

2 2 3 2 2 4

= Жь ^2 = Ж1Ж2, ^3 = Жь ^4 = Ж1Ж2, = Жь

32 4 5 42 24

^>6 = Ж1Ж2, = Ж1Ж2, = Ж1; ^>9 = Ж1Ж2, ^>10 = Ж1Ж2, ...

В качестве пробных функций выбирались те же самые элементы базиса: = фг, ^2 = <Рг- Подставляя данные представления в уравнение (13) и пнтегрпруя по области Б, получаем систему линейных алгебраических уравнений относительно неизвестных коэффициентов разложения а, (¿, г = 1,..., п, решение которой дает искомые коэффициенты разложения функций перемещения.

4.2. Исследование обратной задачи. В качестве дополнительной информации служат данные о п1(ж1 , Ж2), Н2(ж1,Ж2) в некотором наборе точек в области область Б, полученные в результате решения прямой задачи, описанной в п. 4.1. Данная задача является некорректной, так как входная информация о функциях представлена в виде некоторого набора данных, в связи с чем обратная задача решается в два этапа. На первом этапе по обоим наборам значений П1(Ж1,Ж2) и п2(ж1 , Ж2) строится двумерный интерполяционный многочлен Лагранжа; тем самым преодолевается некорректность вычисления производной от функций, заданных таблично. На втором этапе используется проекционный метод, описанный выше; представим искомые функции в виде разложения

п п

А = ^ Сгфг, ^ = ^ (фг, (16)

г=1 г=1

-1 = 1, -2 = Ж1, фз = Ж2, -4 = Ж2, -5 = Ж4, фб = Ж4, ...

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

5. Результаты вычислительных экспериментов

В этом разделе представлены результаты вычислительных экспериментов по решению прямой и обратной задач. Рассмотрена прямоугольная область с параметрами: Б = {ж1 е [0,1], ж2 е [—Ь, Ь]}, где 1 = 100 см, Ь = 10 см, р = 7856 кг/см3. Будем считать, что грань ж1 = 0 жестко защемлена, на грани ж1 = 1 задана касательная нагрузка Р \ _ 1 = 3000 кгс/см2, а грани Ж2 = ±Ь свободны от нагрузок.

5.1. Прямая задача. Для оценки точности полученных результатов решения прямой задачи, проведено сравнение этого решения с решением, полученных с помощью метода конечных элементов. Для оценки точности проанализирована относительная погрешность 6г = где щ — решение задачи проекционным методом, щ — методом конечных элементов.

На рис. 1-3 представлены результаты решения прямой задачи для неоднородного закона изменения А(жьж2) = А0 • (0.3 + + (х)2) и Мжъжг) = Ц-о • (0.3 + (х)2 + (^")2)> где А0 = 0.462 ■ 106 кг/см2, = 0.69 ■ 106 кг/см2. Графически результаты решения П1(Ж1,Ж2), и2(ж1,Ж2) представлены в двух сечениях пластины: при Ж2 = Ь и при Ж1 = I.

Б

де поверхности. На графиках сплошной линией показано решение, полученное методом конечных элементов.

П1 (х1 ,ь)

П2 (XI ,Ь)

П1 (1,Х 2 )

зщ (XI ,х2 )

П2 (1,Х2 )

зщ (XI ,х2 )

Рис. 1. Сравнительные результаты решений прямой задачи проекционным методом и МКЭ в статике (ш = 0). Сверху представлены результаты для 41(х1,Х2), снизу — для 42(х1, Х2).

м1 (х1 ,ь)

м2 (х1 ,ь)

П1 (1,Х2 )

8п, (Х1 ,х2 )

м2 (1,х2 )

Зщ (Х1 ,Х2 )

1П 1ПП

Рис. 2. Сравнительные результаты решений прямой задачи проекционным методом и МКЭ при частоте ш = 3 ■ 2п рад/с.

йщ (XI ,Х2 )

Рис. 3. Сравнительные результаты решений прямой задачи проекционным методом и МКЭ при частоте ш = 11 ■ 2п рад/с.

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

5.2. Обратная задача

Вычислительные эксперименты по реконструкции проводились для разных законов изменения модулей упругости, включая монотонные и немонотонные законы; в данной статье представим лишь некоторые из них. Эксперименты проведены для разных частотных диапазонов до третьего резонанса. Выявлено, что наиболее благоприятным с точки зрения точности реконструкции частотным диапазоном является первый, т. е. до первого резонанса. Для оценки точности восстановления проанализирована относительная погрешность реконструкции 6^ и 6д [11].

На рис. 4-6 приведены результаты по восстановлению неизвестных Х(х\,Х2) и , Х2) та часто те и = 3 • 2п рад/с (первый частотный диапазон). На графиках сплошной линией изображено точное решение, точками — результаты реконструкции; поверхность справа иллюстрирует погрешность идентификации.

,Х2 )

ц,(х1, Ь)

1.5x10

1.4x10®

1.3x10® 6

1.2x10 1.1x10® 1.x 10® 9. х 105

«

•

• •

.1 • •

А(х1, Ь)

, г 1.5x10 Б

1.4x10 , , „п6

1.3x10

1.2x10 6

\ 1.1x10 • 1 ю' • /

1. х 10

Х2)

•

•

/

•

•

Х(1, Х2)

20 40 60 80

• •

\ 1. х 10 „ ,„5- /

\ 9. х 10 5 /

• 8. х 10 •

$п2 (XI ,Х2)

Рис. 4. Результаты реконструкции законов неоднородности Х(хих2) = Ао • (о.З+ + // (./•;. ./••_>) = Аад • (о.З+ +

Ао = 0.462 ■ 106 кг/см2 и [о = 0.69 ■ 106 кг/см2, п = 4.

Анализируя полученные численные результаты, необходимо отметить, что при восстановлении законов неоднородности изменения модулей по некоторым полиномиальным законам, схема работает достаточно эффективно; относительная погрешность не превосходит 0.04% (без учета зашумления входных данных). Для более сложных законов погрешность реконструкции доходит до 9.5%.

Заключение

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

Литература

1. Ватульян А. О. Обратные задачи в механике деформируемого твердого тела.—М.: Физматлит, 2007.-223 с.

2. Ватульян А. О. К теории обратных задач в линейной механике деформируемого тела // ПММ.— 2010.—Т. 74, № 6.-С. 909-916.

3. Bonnet М., Constantinescu A. Inverse problems in elasticity // Inverse Probl.—2005.—№ 21.—P. 1-50.

4. Бочарова О. В., Ватульян А. О. О реконструкции плотности и модуля Юнга для неоднородного стержня // Акустический журн.—2009.—Т. 55, № 3.—С. 275-282.

5. Санчес-Паленсия Э. Неоднородные среды и теория колебаний.—М.: Мир, 1984.—472 с.

6. Тихонов А. Н., Арсенин В. Я. Методы решения некорректных задач.—М.: Наука, 1986.—287 с.

7. Филиппов А. 17. Колебания деформируемых систем.—М.: Машиностроение, 1970.—736 с.

8. Гюнтер Н. М. Интегрирование уравнений первого порядка в частных производных Ленинград: Гостехиздат, 1934.—359 с.

9. Dudarev V. V., Nedin R. D., Vatulyan A. O. Nondestructive identification of inhomogeneous residual stress state in deformable bodies on the basis of the acoustic sounding method // Advanced Materials Research.-2014.-Vol. 996.-P. 409-414.

10. Nedin R. D., Vatulyan A. O. Inverse Problem of Non-homogeneous Residual Stress Identification in Thin Plates // Int. J. Solids Struct—2013—№ 50.-P. 2107-2114.

11. Ватульян А. О. Гукасян Л. С. О задаче Коши для уравнения в частных производных первого порядка и ее приложениях в теории обратных задач // Вестн. ДГТУ.—2012.—Т. 68, № 7.—С. 1120.

Статья поступила 15 сентября 2015 г.

Ватульян Александр Ованесович Южный федеральный университет, заведующий кафедрой теории упругости РОССИЯ, 344090, Ростов-на-Дону, ул. Мильчакова, 8 а Южный математический институт ВНЦ РАН, заведующий отделом дифференциальных уравнений РОССИЯ, 362027, Владикавказ, ул. Маркуса, 22 E-mail: vatulyan@aaanet.ru

Гукасян Лусинэ Суреновна

Донского государственного технического университета, старший преподаватель кафедры прикладной математики РОССИЯ, 344010, г. Ростов-на-Дону, пл. Гагарина, 1 E-mail: luska-90@list.ru

Недин Ростислав Дмитриевич

Южный математический институт ВНЦ РАН,

младший научный сотрудник

РОССИЯ, 362027, Владикавказ, ул. Маркуса, 22

E-mail: rdn90@bk.ru

ON THE CAUCHY PROBLEM IN THE THEORY OF COEFFICIENT INVERSE PROBLEMS FOR ELASTIC BODIES

Vatulyan A. O., Gukasyan L. S., Nedin R. D.

The inverse problems on the identification of inhomogeneous material properties of an elastic plate is considered. The condition of uniqueness of the inverse problems statements are analyzed. Direct problem on finding displacements to formulate the additional input data of the inverse problem is investigated; the accuracy of the direct problem solution is estimated by means of comparison with finite-element computation. A scheme of the inverse problem solving is proposed based on the application of the weak statement of the initial boundary problem and the projection method. A series of computation experiments on a reconstruction of various types of inhomogeneity laws of the Lame' coefficients is performed.

Key words: coefficient inverse problem, the Lame' coefficients, weak statement, projection method, finite element method, biharmonic polynomials.