Граничные интегральные уравнения первого рода с гладкими ядрами в динамических задачах ортотропной теории упругости

Рассмотрим установившиеся гармонические колебания ортотропного тела в области Ω ⊂ R ² с кусочно-гладкой границей ∂Ω . Оси ортотропии совпадают с осями x ₁ ,x ₂ . Временна´я зависимость имеет вид: e ^{- iωt} , где ω — циклическая частота колебаний, t — время, i — мнимая единица (г ² = -1 ). Запишем закон Гука при плоском напряжённом состоянии [12] :

C 22

где e ij = ( u i,j + U ji ) / 2 , и уравнения движения:

f C11 u 1 , 11 + ^C 66 u 1 , 22 + (C 12 + C 66 ^{) u} 2 , 12 +pw ² u i =^{0 ,} [ (^C 12 + C 66 ^{) u} 1 , 12 + ^C 66 u 2 , 11 +C 22 u 2 , 22 + pW 2 u 2 =0.

Введём в рассмотрение вектор напряжений с компонентами на границе t j = a ij n i , где a ij = C ij _kl u _k ,l , и сформулируем общую смешанную граничную задачу:

дп =S 1 US 2 ,     u j |_Si = f j ,

^{t j |} s ₂ = g j .

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

• 3.    Вывод альтернативных граничных интегральных уравнений — АГИУ

  • 3.1.    Спектральная форма

Введём для компонент вектора перемещений преобразование Фурье по пространственной переменной x ∈ Ω :

U

j ^{( a )} = /

u j

( x ) e i ^a'x dA,

a = (a 1 ,a ₂ ) GR ² .

Ω

Преобразование Фурье уравнений (2) с учётом формулы Грина приводит к матричному уравнению относительно трансформант

M - ( a ) u ( a ) = F ( a ),                                             (5)

в котором характеристическая матрица симметрична и имеет вид:

ρω ² - C 11 α ² ₁ - C 66 α ² ₂ \ ^{—( C} 12 +C 66 ^{) a} 1 ^a 2

^{— ( C} 12 + C 66 ^{) a} 1 ^a 2 ρω ² - C 66 α ² ₁ - C 22 α ² ₂

а компоненты вектора правой части (5) выражаются через граничные операторы. Введём обозначение:

Bjk (a,n) iapCpjklnl, где суммирование выполняется по повторяющимся индексам p,l G {1,2}. Тогда

F j ( a ) = I [ B jk ( a , n ) u k - t j ] e i ^{a x} dS.

∂Ω

Здесь суммирование ведётся по повторяющемуся индексу k G {1,2} , индекс j G {1,2} является свободным.

Подстановка известных граничных условий в (7) даёт явную связь с условиями (3) . В смешанной постановке неизвестны t j на S ₁ и u j на S ₂ :

F j ( a )= / [B jk ( a , n ) f k -t j ] e i ^{a x} dS + / [ B jk ( a , n ) u k - g j ] e i ^{a x} dS.

S 1

S 2

Случай S 1 = дП соответствует первой граничной задаче теории упругости (неизвестны t j ), случай S 2 = дП — второй краевой задаче (неизвестны u j ). Форма АГИУ при этом сохраняется, изменяется лишь состав неизвестных.

• 3.2.    Характеристическое множество

Компоненты перемещений будем рассматривать как достаточно регулярные функции, отличные от нуля лишь в ограниченной области Ω и продолженные нулём вне её. По теореме Пэли–Винера [13] преобразования Фурье (4) таких функций с вещественных значений α ∈ R ² распространяются на комплексные α ∈ C ² и оказываются целыми функциями переменных α 1 , α 2 . Они аналитичны всюду (наподобие многочленов и экспоненты), а значит, не имеют ни полюсов, ни иных особых точек. Таким образом, вектор трансформант и ( а ) свободен от особенностей при любых α ∈ C ² . Именно это свойство лежит в основе дальнейшего вывода. Обратную матрицу удобно представить через присоединённую adjM _w . Тогда M -1 = adjM _w /detM ^ , и формальное обращение системы (5) записывается как u adjM ^ ( а ) F ( а )/detM _ш ( а ) . У правой части возникают полюса в точках, где del М ^ ( а ) = 0 . Устранить это противоречие для и можно лишь при обращении в нуль числителя этой дроби — вектора adjM_w ( a ) F ( a ) — на этом множестве.

Введём в рассмотрение характеристическое множество:

К _ш = { а е C 2 :detM _w ( а ) = 0}.

Для ортотропного материала его структура определяется уравнением:

^detM _w = р 2 Ш 4 р^ 2 [ ⁽С 11 + C 66 ^{) a} 1 + ^{( C} 22 + ^C 66 ^{) a} 2 ] +D 4 ,                      ⁽⁹⁾

где

D 4 = C ii C 66 a 4 + (C_UC 22 - C 22 -2С 12 С бб )а 2 а 2 +С 22 С бб а ⁴ .                    (10)

• 3.3.    Спектральные условия

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

С учётом введённого в подразделе 3.2 представления обратной матрицы через присоединённую условие разрешимости системы (5) с матрицей (6) на характеристическом множестве K ω принимает вид:

adjM _w ( a ) F ( a )= 0 ,     а ■ K - ,                                   (11)

где для скалярной формы записи элементы матрицы (6) обозначимкак M 11 ,M 22 ,M 12 . Тогда условие разрешимости становится следующим:

M 22 F 1 -M 12 F 2 = 0,    -M 12 F 1 + M 11 F 2 = 0,     а еК _ш .

В регулярной точке характеристического множества ранг матрицы М ^ ( а ) равен единице. Для матрицы второго порядка это означает, что ранг присоединённой матрицы adjM _w ( а ) также равен единице, поэтому две скалярные строки в (12) линейно зависимы, а условие (11) содержит только одно независимое скалярное уравнение.

Параметризуем множество К _ш набором компонент a ^(s) (т) , где т — параметр вдоль соответствующей ветви характеристического множества. Требование выполнения условия (11) на каждой компоненте приводит к системе АГИУ 1-го рода с гладкими ядрами:

adjM _w ( a ^(s) (т) ) F ( а ^(s) (т) ) = 0 (s = 1,2). (13)

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

В смешанной постановке учёт в скалярных условиях (12) соотношения (8) даёт систему АГИУ относительно неизвестных t j на S 1 и u j на S 2 :

X (adj M _w ) mj tlx ( x^ i^ ' x

j =1

dS — У B jk ( а , п )и _к ( x^i^dS

S 2

= E(adjM )mj / j=1               S1

B jkM f k ( х ) е^м' ^х dS — ^j ( x^dS

S 2

( m =1,2),      а = а ^(s) (т).

Система АГИУ (14) допускает различные варианты дискретизации. При этом термин «дискретизация» употребляется в двух различных смыслах:

- дискретизация границы — разбиение ∂ Ω на участки и аппроксимация граничных функций конечным числом коэффициентов;

• - дискретизация характеристического множества — выбор конечного набора точек на его компонентах.

• 3.4.    Параметризация

В настоящей работе идёт речь о дискретизации характеристического множества. Используется коллокационный подход, при котором операторные уравнения (14) удовлетворяются в выбранных точках α ∈ K ω . Альтернативой служат проекционные методы, в частности метод Галёркина.

Приведём структуру характеристического множества для некоторых упругих материалов. Для изображения характеристического множества введём безразмерные спектральные координаты ^ j = R 0 a j (j = 1,2 ), где R 0 — характерный размер области, α j — компо нент ы спектрального параметра преобразования Фурье (4) . Частоту зададим в безразмерном виде: А _ш = wR 0 у/p/E 2 .

На рисунке 1 представлены вещественные ветви К _ш характеристического множества при А _ш = 1.2 . Внешняя ветвь является выпуклой, а внутренняя в зависимости от параметров материала может терять выпуклость. Кривые построены для материалов, рассматриваемых в главных осях ортотропии. Использованы безразмерные параметры x E = E 1 /E 2 , X g = G 12 /E 2 , v = v ₁₂ . Здесь E 1 и E2 — модули Юнга вдоль главных осей ортотропии x 1 и x 2 ; G 12 — модуль сдвига в плоскости ( x 1 ,x 2 ); ν 12 — коэффициент Пуассона, характеризующий поперечную деформацию вдоль оси x 2 при одноосном нагружении в направлении x 1 . В изотропном случае обе ветви являются концентрическими окружностями (Рис. 1a ).

Рис. 1. Вещественные ветви характеристического множества (внешняя ветвь – сплошная линия, внутренняя ветвь – пунктирная линия) для материалов с разными свойствами: изотропный ( χ E = 1 , χ G = 0 . 385 , ν = 0 . 3 ) ( а ); ортотропный ( χ E =3 , χ G = 0 . 5 , ν = 0 . 25 ) ( б ); древесина ( χ E =20 , χ G = 1 . 2 , ν = 0 . 35 ) ( в )

Характеристическое множество допускает разные параметризации. Далее используется представление в полярных координатах. При а1 = a(s) cos0, а2 = a(s)sin0 уравнение (9) определяет обе ветви как a«W =           ±              ,    (s = 1,2), где

a(0) = C11C66 cos40+(C11C22-C^2 - 2Ci2C66)cos2 0sin20+C22C66sin40,(16)

Ф) = (C11 + C66)cos20 + (C22+C66)sin20.

Заметим, что альтернативой может служить декартова параметризация. При фиксировании одного из параметров ( a 1 = £ или a ₂ = n ) корни биквадратного уравнения (9) задают соответственно точки или одной группы ветвей, или другой. Такая параметризация больше подходит для областей с прямолинейными участками границы. Граничные интегралы содержат множители e ^{iα 1 x 1} и e ^{iα 2 x 2} и допускают явное покоординатное интегрирование. Независимо от вида параметризации возможен различный выбор точек коллокации на K ω , в которых выполняются условия АГИУ (13) .

• 3.5.    Статический предел

Приведённый статический переход относится к ортотропным материалам с различными корнями характеристического биквадратного уравнения. При ш = 0 аналогичные рассуждения приводят к статическому варианту АГИУ. При этом характеристическое уравнение (9) сводится к однородному уравнению 4-й степени:

D4 (а) = С11Сбба1+Ва1а2+С22Сбба2 = 0, B = C11C22-C222 - 2C12C66, которое определяет статическое характеристическое множество ненулевых комплексных спектральных векторов а = (а1,а2):

К о = { а е C ² \{ 0 }: D 4 ( а ) = 0}.

Здесь C ² — пространство упорядоченных пар комплексных чисел, а 0 = (0,0) — нулевой вектор. Таким образом, условие а е C ² \{ 0 } означает, что компоненты а 1 ,а ₂ комплексные и не равны нулю одновременно.

Положительная определённость упругой энергии исключает ненулевые вещественные решения уравнения D 4 ( a ) = 0 , то есть К _о ПК ² = 0 . Покажем это, положив в выражении (6) ш = 0 . Согласно (9) , detM ₀ ( a ) = D 4 ( a ) . Для проверки положительной определённости матрицы M 0 ( а ) возьмём произвольный ненулевой вещественный спектральный вектор а е R ² и произвольный ненулевой постоянный вектор а = (а 1 ,а ₂ ) е R ² . Рассмотрим вспомогательное линейное поле перемещений u * ( x ) = а ( а • x ) . По определению малых деформаций, приведённому после выражений (1) , этому полю отвечает постоянный тензор:

*     u i,,j + u j,i    a i ^a j + a j ^a i

( i,j = 1,2).

^e j =---2---=----2----

И это есть тензор деформаций вспомогательного поля, а не его преобразование Фурье. Вектор a задаёт направление и величину перемещений этого поля. Плотность упругой энергии W (е * ) = C ijki e j e kl / 2 связана с характеристической матрицей равенством:

а ^T [-M о ( a )] а = C _ijk i e *j е^ =2W (е * ) > 0.

Здесь суммирование выполняется по повторяющимся индексам от 1 до 2. При ненулевых вещественных a и α тензор ε ^∗ ненулевой, поэтому строгое неравенство вытекает из положительной определённости упругой энергии. Следовательно, матрица М ₀ ( а ) положительно определена, а её определитель D 4 ( а ) положителен при любом ненулевом вещественном а (см. работу [3] ). Подстановка а ₂ = ^а 1 приводит к биквадратному уравнению:

C22C66M4 + BM2 + C11C66 = 0, четыре корня ^r которого задают комплексные характеристические прямые a2 = ^rа1. В статике коллокационные точки берутся на комплексных прямых

а(г)(к) = (к,^г к),     к е R\{0}     (r =1,...,4), где µr — корни статического уравнения для µ, приведённого выше. При декартовой параметризации к ним можно добавлять симметричное семейство. Таким образом, в статическом пределе сохраняются структура спектральных условий и набор граничных неизвестных.

Покажем, откуда вытекают первые интегральные следствия этой предельной формы. На характеристической прямой а(г)(к) = (к,цгк) устремим к ч 0 и разложим экспоненту в (7): ei^x = 1 + iа • x + О(к2). При этом операторы Bjk (а,п) линейны по а и исчезают при а ч 0. Приравнивание нулю коэффициентов при младших степенях κ в спектральном условии даёт: в нулевом порядке — равенство нулю главного вектора поверхностных усилий, а в первом порядке — главного момента, то есть классические условия равновесия сил и моментов в плоской задаче:

jtdS = 0 ,

∂Ω

У (x 1 t ₂ — x ₂ t 1 ) dS = 0.

∂Ω

• 4.    Численная реализация и регуляризация

  • 4.1.    Дискретизация

Схема решения базируется на граничной дискретизации. Кусочно-гладкая граница области разбивается на элементы:

P дП = J Гр,    Гр = {x = Xp(Z ):0 < Z < 1}    (Р =1,...,P), p=1

где P — число геометрически гладких участков границы, а ζ — локальный параметр на Γ p . Для каждого участка определена внешняя единичная нормаль n = n _p (Z) , входящая в граничные операторы (7) , (8) .

Состав граничных неизвестных задаётся постановкой (3): в смешанной краевой задаче теории упругости на участках S1 неизвестны компоненты вектора напряжений tj, а на участках S2 — компоненты вектора перемещений uj (j = 1,2). Эти неизвестные удобно записать единообразно: обозначим через qpk(Z) (k = 1,...,mp) скалярные неизвестные граничные функции на участке Гр — компоненты вектора напряжений при Гр С S1 и компоненты вектораперемещений при Гр С S2. Независимыми считаются функции, остающиеся после учёта связей, налагаемых симметрией задачи и условиями согласования в угловых точках (подробнее показано в подразделе 5.2). Тогда число P скалярных блоков равно m =

m p .

В частности, при двух неизвестных компонентах на каждом участке m = 2 P .

p =i

При сквозной нумерации блоков l = 1,... ,m через N l обозначим число базисных функций в l -м блоке. Переменные в каждом блоке аппроксимируем линейной комбинацией:

Nl                                                           m qi(Z) ~ qNKZ ) = '^JwlnФln(Z),     w = (wii,...,wiN1 ;...;Wm1,...,WmNm )T,     N = ^Nl.

n =1                                                                 l =1

Здесь N — полное число коэффициентов всех неизвестных, а w — вектор коэффициентов. В расчётах в качестве ϕ ln используем полиномы Чебышёва 1-го рода, заданные на всём участке Γ p . В отличие от функций формы МКЭ такой базис глобален: его функции не имеют локального носителя, а коэффициенты w ln не являются узловыми значениями искомой функции. Значение последней в любой точке восстанавливается всей суммой разложения. Подстановка такого представления в (8) даёт равенство:

N

F ( a ; w ) = F ⁽⁰⁾ ( a ) + ^ W n G n ( a ), (19) n =1

где F ⁽⁰⁾ содержит вклад заданных функций f j , g j , а G _n — вклад одной базисной функции неизвестного перемещения или напряжения.

Переход к конечномерной системе осуществляется выбором точек α q на ветвях K ω . В статическом случае используются точки на К ₀ , заданные комплексными прямыми. В каждой точке коллокации две строки условия (12) линейно зависимы (см. подраздел 3.3), поэтому в дискретной системе достаточно оставить одну из них, а именно берётся строка r * ( q ) с большей евклидовой нормой (из двух пропорциональных строк она вычисляется с меньшей относительной погрешностью):

A C n = [adjM _w ( a _g ) G n ( a ₉ )] r _{( q )} , f C = - | adjM _w ( a ₉ ) F ⁽⁰⁾ ( a ₉ ) ] (20) L r * ⁽ 9)

Верхний индекс C указывает, что определённые в (20) элементы A _q ^C _n и компоненты правой части f _q ^C являются комплексными числами. Составленные из них матрицу и вектор обозначим соответственно через A C и f C . При декартовой параметризации вводятся сетки {^ i } N ^£ ₁ и {n i }^ . Каждому фиксированному значению ^ i соответствуют четыре корня a j)^, ) (j = 1,...,4 ) характеристического уравнения. Аналогично, при фиксированном n i находятся четыре корня a 1 ^j) (n i ) . При учёте всех корней получается число Q = 4 N +4 N _n комплексных точек коллокации. Для вещественных коэффициентов w _n комплексная система расщепляется на вещественную и мнимую части:

Aw = f , A = f^ A CY f = (^ C ) , A e R ^2QxN . (21) im A ши

Здесь Q — число комплексных коллокационных точек. После разделения на вещественную и мнимую части получается 2Q вещественных уравнений. В результате имеем переопределённую систему с матрицей A e R ^{2 QxN} , вектором неизвестных w e R N и правой частью f e R ^2Q . В общем случае количество базисных функций N l может зависеть от блока. Так, при равномерной по блокам аппроксимации полагаем N l = N c . Тогда при двух независимых скалярных функциях на каждом участке будем иметь: N = 2 PN c . Заметим, что эта формула относится к системе до возможной симметрийной редукции и наложения дополнительных связей между коэффициентами.

Воспользуемся схемой избыточной коллокации, согласно которой число вещественных уравнений выбирается существенно большим числа неизвестных ( 2Q >  N ). Поскольку число уравнений значительно превышает число неизвестных, прямоугольная матрица A имеет ранг rank A <  N . Подчеркнём, что выбор одной строки в (20) переопределённости не устраняет, а лишь исключает дублирующее (линейно зависимое) уравнение в каждой точке коллокации, тогда как система (21) остаётся переопределённой и решается в смысле наименьших квадратов. Необходимость регуляризации, о которой пойдёт речь ниже, вызвана не переопределённостью, а быстрым убыванием сингулярных чисел матрицы A .

• 4.2.    Критерий усечённого сингулярного разложения (УСР) и параметры регуляризации

Гладкость ядер, дающая формулировке АГИУ её аналитические преимущества, имеет вычислительную цену: оператор с гладким ядром сглаживает любые граничные данные, поэтому обратная к нему задача — восстановление граничных функций — чувствительна к малейшим возмущениям. На дискретном уровне неустойчивость проявляется в быстром убывании сингулярных чисел матрицы A. Для подавления неустойчивости применяются усечение сингулярного разложения и регуляризация А.Н. Тихонова [14]. При УСР вклады, отвечающие малым сингулярным числам, отбрасываются; регуляризация Тихонова ограничивает решение посредством стабилизирующего функционала.

В классическом методе А.Н. Тихонова для системы Aw = f минимизируется стабилизирующий функционал:

llAw-fIIL2 + A|Lw|L2 ^min,    A>0, где L — оператор регуляризации, выбор которого задаёт априорные требования к искомому решению. В простейшем случае при L = I ограничивается норма решения: при учёте гладкости вместо единичного оператора вводится стабилизирующий оператор более высокого порядка, соответствующий норме производных. Его дискретное представление может строиться на основе конечных разностей либо, как в используемых ниже базисных аппроксимациях, путём непосредственного дифференцирования базисных функций.

С учётом блочной структуры неизвестных классическая схема А.Н. Тихонова принимает вид:

m

A i > 0.

|Aw — f ^2 +XMLlwl|lL2 ^ min, l=1

Здесь используется введённый выше блочный вектор неизвестных w =( w 1 ,..., w _m ) ; L l — оператор регуляризации в l -м блоке и λ _l — соответствующий параметр регуляризации. Состав блоков, числа базисных функций и операторы L l выбираются с учётом параметризации границы, типа неизвестных граничных функций и предполагаемой гладкости искомого решения.

Параметры регуляризации выберем по внутренним характеристикам дискретной системы (21) , что соответствует стандартной практике решения дискретных некорректных задач [14, 15] . Пусть

A = USV T ,     S = diag(a_b...,a N _a ),    a i >  ••• >>  0

есть компактное сингулярное разложение матрицы. Здесь столбцы U k и V k матриц U и V — соответственно левые и правые сингулярные векторы, а Σ — диагональная матрица сингулярных чисел σ k . Сингулярные числа вычисляются как неотрицательные квадратные корни собственных значений матрицы A ^T A и определены для любой матрицы. Компоненты решения, отвечающие малым σ k , в наибольшей степени усиливают погрешность правой части. На отбрасывании таких компонент и основано усечение, давшее название методу. Число N = min(2Q,N) задаёт количество сингулярных чисел дискретной матрицы. Для используемых далее систем обычно N _a = N . В УСР удерживаются только первые r сингулярных чисел, где 1 <  r <  N [15] , а приближённое решение имеет вид:

w r

r

=X k=1

( U k , f ) σ k

V k .

Рабочий ранг r определяется по относительному спектру σ k /σ 1 после отсечения малых сингулярных чисел. На практике используется относительный порог ε _svd и, при наличии выраженного спектрального разрыва, выбирается такой r , что

σ r σ r +1

⩾ ε svd , < ε svd .

σ 1 σ 1

Параметры блочной регуляризации λ l выбираются после фиксации масштабирования системы и состава блоков неизвестных. Они подавляют неустойчивые высокочастотные компоненты по принципу регуляризации А.Н. Тихонова [14] . В настоящей работе для автоматического выбора λ l используется критерий обобщённой перекрёстной проверки — Generalized Cross-Validation (GCV) [16] . При обычной перекрёстной проверке поочерёдно исключают одно уравнение системы (21) , находят решение по оставшимся уравнениям и оценивают невязку в исключённом уравнении. В GCV индивидуальные поправки, учитывающие влияние каждого уравнения на решение, заменяются общей поправкой через след матрицы влияния. Поэтому критерий вычисляется непосредственно по всей системе (21) , без поочерёдного удаления её строк.

Обозначим через A + G R N ^{x 2Q} оператор регуляризованного решения всей системы (21) : w _A = A + f . Здесь λ обозначает параметр регуляризации либо совокупность блочных параметров λ l при фиксированном рабочем ранге УСР. Матрица Н _А = AA + G R ^{2Q x 2Q} отображает правую часть f в вектор восстановленных значений Aw _a и называется матрицей влияния. Критерий GCV, с точностью до постоянного множителя, не влияющего на положение минимума, имеет вид [16] :

G(A)= ^{| A} w ^{A f 1} , ,     H_A = AA + ,     I GR ^{2Q x 2Q} .

[ tr( I - Ha )]^{2                 A}

Во всех частях этой формулы присутствует одна и та же матрица A системы (21) , содержащая все её 2Q строк. Оператор A ⁺ построен также без удаления проверяемого уравнения. Выбор одной из двух линейно зависимых строк условия (12) уже выполнен при формировании матрицы (см. (20) ) ине связан с процедурой перекрёстной проверки. Минимум G определяет параметры регуляризации без априорной оценки уровня погрешности правой части.

Дискретизация контролируется по устойчивости результата: при умеренном изменении Q , распределения точек на K ω и числа базисных функций эффективный ранг, структура сингулярного спектра и восстановленные граничные функции должны сохраняться. В некоторых задачах имеются особенности граничных полей напряжений в угловых точках, которые могут быть установлены из рассмотрения однородной граничной задачи теории упругости для ортотропного клина. К особым точкам относятся угловые точки кусочно-гладкой границы, точки смены типа граничных условий и разрывы заданных граничных данных. Однако локальные особенности граничных функций не изменяют вид спектрального условия (11) , но должны учитываться при выборе базиса. Для изотропной теории упругости исходной является классическая асимптотика Вильямса [11] ; в анизотропных и составных клиньях используются её обобщения [3, 17, 18] .

В окрестности особой точки x c введём локальную полярную систему координат ( r,0 ) с полюсом x c в вершине клина, расстоянием до вершины r и углом θ , отсчитываемым от одной из сторон клина внутрь его области. Если в этих координатах у вектора перемещений главный член представляется в виде:

u(r,0) ^rvФ(0),     r = |x—xc|, то напряжения имеют порядок rv-1. Энергетически допустимыми являются решения с Rev > 0. При 0 < Rev < 1 напряжения могут быть неограниченными, однако принадлежность решения к энергетическому классу сохраняется. Эта особенность должна учитываться в локальном базисе на элементе в окрестности особой точки.

В дискретной схеме около каждой особой точки выделим локальный участок длиной δ . Для неизвестной граничной функции q используем представление:

q(r) = q^Tes ( r )+ r^Y ^d m ^ _m ( r/S ) ,

m где qreg — регулярная часть, ^m — гладкие функции нормированной локальной координаты, а y = Rev — 1 — показатель особенности для компонент вектора напряжений на границе. Для спектральных задач рассмотрим однородную граничную постановку. Собственную частоту определим как значение параметра, при котором существует ненулевое поле перемещений, удовлетворяющее уравнениям движения в области и однородным условиям на границе. Для численного поиска воспользуемся введённой в подразделе 3.4 безразмерной частотой:

Л . = uR ₀ ^Р/Ё ₂ ,

где R ₀ — характерный размер области. Резонансные значения Л _ш найдём по спектральному индикатору, основанному на конечномерной аппроксимации однородной задачи методом частных решений (МЧР). По своей конструкции такой подход близок к методам Ритца и Треффтца: искомое поле ищется в конечномерном пространстве пробных решений. Использование полных несингулярных функций Треффтца для гармонических и статических задач плоской упругости развивалось в работах [19, 20] . Близкая по конструкции схема с экспоненциальными базисными функциями, параметры которых определяются характеристическим уравнением оператора, применялась к статическим и гармоническим задачам упругости [21] . Однако, в отличие от классической схемы Ритца, пробные поля здесь точно удовлетворяют уравнениям движения в Ω , в сравнении же с произвольным базисом Треффтца их параметры выбираются на характеристическом множестве K ω , что согласуется с АГИУ-постановкой [5 –10] . Стабилизация индикатора соответствует современному варианту МЧР [22, 23] .

Для фиксированного λ ω базис пробного пространства составляют экспоненциальные решения (2) . Пусть a m ) (Л _ш ) = (a l m (Л _ш ) ,a (, m (Л _ш )) . Теперь при построении пробного пространства, в отличие от сквозной нумерации точек коллокации в (20) , используем двойную нумерацию: s — номер ветви характеристического множества, m — номер выбранной точки на ней, а a m ) (Л_ш) — собственный вектор характеристической матрицы (6) . Базисные поля перемещений имеют вид:

v _m s ( x 1 ,X 2 ; Л _ш )= a m )^ )exp{i^a l m (Л . )x i + a , m (Л . )x ,] }.

Здесь vms — двухкомпонентное векторное поле перемещений. Общее число выбранных базисных полей обозначим через N . Их линейная комбинация со скалярными коэффициентами cms записывается как vc (x1 ,x2 ;Л. )  ^^cmsvms (x 1 ,x2 ;Лш ),     c  {cms }.

m,s

Вектор напряжений для пробного поля не вводится как независимый параметр, а вычисляется непосредственно по закону Гука (1) и соотношению t j = n i , приведённому после (2) . Зависимость пробных полей от констант материала учитывается через характеристическую матрицу (6) , поэтому и точки a m ) (Л _ш ) , и векторы a m ) (Л _ш ) определяются упругими константами материала.

Собственная частота рассчитывается по минимуму стабилизированной невязки однородных граничных условий:

\      ^R ^A “ v c ^2 ( d^ )

n N (А _ш ) = mm—-— -------,     Ф ( Х _Ш ) = -log ₁₀ n N (А _ш ),                    (24)

c = ⁰     ll v c || L 2( n )

ll ^v c ^| L 2( n ) = У ^{| v} c ( ^{x ) | 2 d} A.

Ω

Минимум в (24) берётся по коэффициентам c , для которых | v _c | L ₂ ( ^ ) > 0 . При линейно независимом базисе это равносильно неравенству: c = 0 . Оператор R x ^ задаёт невязку однородных граничных условий для пробного поля v _c . Например, для свободной границы она совпадает с компонентами вектора напряжений на границе; на гладком участке при единичной нормали n и единичном касательном векторе s можно записать: R x ^ v _c = (p _n ,p _s ) . Здесь p j = a ij n i , p _n = p • n , p s = p • s , где j — компоненты тензора напряжений Коши для пробного поля, n и s — единичные нормаль и касательный вектор к границе. Нормировка на знаменатель делает отношение не зависящим от масштаба c . Кроме того, она исключает ложные минимумы, связанные с почти тривиальными линейными комбинациями, и соответствует стабилизированному варианту МЧР [23] .

Размер пробного пространства выбирается по стабилизации спектрального индикатора: при последовательном расширении пробного пространства положения локальных минимумов η N и форма соответствующих пиков Φ на резонансах должны изменяться слабо.

• 4.3.    Симметрия и классы собственных частот

• 5.    Вычислительные эксперименты

• 5.1.    Эллиптическая область с заданным вектором напряжений

Если геометрия области, упругие коэффициенты и граничные условия инвариантны относительно отражений x 1 н —x 1 и x ₂ H —x ₂ , то дифференциальный оператор коммутирует с действием группы Z ₂ х Z ₂ . В этом случае пространство допустимых полей раскладывается в прямую сумму четырёх инвариантных подпространств, соответствующих неприводимым представлениям этой группы [24, 25] .

Класс симметрии задаётся парой знаков (х 1 ,Х 2 ) €{±1} ² : скалярное поле класса (х 1 ,Х 2 ) умножается на х 1 при отражении x 1 н — x 1 и на х ₂ при отражении х ₂ н —х ₂ . Компоненты перемещений не скаляры. При отражении оси направленная вдоль неё компонента тоже меняет знак, поэтому чётности компонент имеют вид: u 1 : (—х 1 ,Х 2 ) , u ₂ : (х 1 ,-Х 2 ) . В дальнейшем введём обозначения: A =(+,+) , B = (—,+) , C =(—,—) , D = (+,—) . Так, в классе A компонента u 1 нечётна по x 1 и чётна по x ₂ , а u ₂ — чётна по x 1 и нечётна по x ₂ . Поиск собственных форм может проводиться независимо в каждом из этих подпространств, а их объединение даёт полный набор симметрийных типов.

Далее рассмотрены решения нескольких задач, иллюстрирующих описанный выше подход. Необходимые для валидации эталонные расчёты выполнены при помощи МКЭ в пакетах FEniCS (с реализацией на сетках из треугольных элементов) и scikit-fem (с элементами типа P ₂ ), а также получен МКЭ-эталон реализации задачи в спектральной постановке. В численных примерах выбранные значения рабочего ранга УСР и локальных параметров регуляризации проверены при вариации параметров дискретизации. Для спектральных задач дополнительно проконтролирована устойчивость положения резонансных значений.

Задача решалась в области П , ограниченной эллипсом Г = д^ :

Г =/(x i ,X 2 ): x 1 + x 2 = 1V    a = 0.125 м , b = 0.080 м .

a ² b ²

Граница параметризовалась как х (ф) = (acosф,bsinф) , 0 <  ф <  2п . На верхней и нижней дугах поверхности, которым отвечают углы 30 ° , нагрузка задавалась в виде нормальной и касательной компонент:

tn(ф) = —Po{q(ф —90°) + q(ф—270°)},    tT (ф)=0,    Ро = 1 МПа, где

q ⁽^ =

^{| s |} < 1,

|s| >  1,

ψ

s =-- 15°

.

При этом Р _о — масштаб нагрузки. Выбранная гладкая функция имела значение q(0) = 2 , а максимальный модуль нормального поверхностного напряжения на каждой нагруженной дуге равнялся 2Р ₀ .

Главные оси ортотропии совпадали с осями эллипса. Исследовался материал — однонаправленный углепластик (близкие справочные значения для системы T300/5208 [26]), обладающий ярко выраженными ортотропными свойствами со следующими упругими константами: E1 = 150 ГПа, E2 = 10 ГПа, G12 = 5 ГПа, ν12 = 0.30, ν21 = 0.02. Этим значениям соответствовали коэффициенты матрицы жёсткости:

/ 150.91 3.02 0 \

C = 3.02 10.06 0 ГПа , ρ=1600 кг / м ³ . 0 0 5.00

Частота ω = 75 00 рад/с отвечала λ ω =ωR 0 ρ/E 2 = 0.3 , где для эллиптической области принимался характерный размер R 0 = ab= 0.1 м.

Неизвестными являлись полярные компоненты перемещения u r , u φ ; масштаб перемещений считался равным U _ref = P ₀ R ₀ /E ₂ = 10 ^{- 5} м. При построении коллокационной системы использовалось по 40 значений каждого координатного параметра характеристического множества: N ξ = N η = 40 . Для каждого фиксированного значения параметра учитывалось четыре корня характеристического уравнения, поэтому число комплексных точек составляло Q = 4N ξ + 4N _η = 320 , что после разделения на действительную и мнимую части дало 640 вещественных уравнений. Число коэффициентов разложения выбиралось в соответствии с критерием устойчивости из раздела 4 (дальнейшее расширение базиса не изменяло восстановленные граничные функции в пределах принятой точности и не увеличивало информативную часть сингулярного спектра). В итоговом расчёте использовалось число коэффициентов N = 18 : N u _r = 10 для u r и N u _φ = 8 для u φ , так что N = N u _r +N u _φ . Для матрицы размером 640 х 18 число сингулярных чисел равнялось min(640,18) = 18 . Устойчивость выбранной конфигурации контролировалась при варьировании N ξ = N _η и числа коэффициентов в обоих блоках; рабочий ранг r = 14 соответствовал устойчивому отделению первых сингулярных значений от области малых сингулярных чисел. Получающаяся деформация являлась симметричной: нагруженные верхняя и нижняя дуги смещались внутрь, а боковые участки — наружу на меньшую величину.

Построенное решение сравнивалось с реализацией посредством МКЭ с локальным сгущением в окрестности нагруженных дуг (сетка состояла из 42 720 элементов и имела 172 074 степени свободы). Погрешность АГИУ составила: 4.6% для u r (по всей границе) и 2.8% для u φ (по всей границе); на дугах — участках приложения симметричной нагрузки — ошибка u φ возрастала до 3.6%, что связано с локализованным характером заданного поверхностного напряжения.

Для контроля устойчивости дискретизации дополнительно анализировался спектр нормированной матрицы системы. Спектральные характеристики служили внутренней диагностикой дискретной системы: по ним выбирался эффективный ранг УСР и проверялось отделение информативных сингулярных компонент от области малых сингулярных чисел. Для компонент u r и u φ использовалась блочная регуляризация с оператором второй производной L = D ² . Аппроксимация выполнялась полиномами Чебышёва 1-го рода T n (τ) , где локальный параметр участка линейно отображается на τ е [-1,1] . Два параметра регуляризации, λ _{u r} и λ _{u φ} , определялись по минимуму GCV-функционала (см. подраздел 4.2), обобщённого на случай двух параметров: G(λ u _r ,λ u _φ ) минимизировался по двумерной сетке значений для блочной регуляризованной системы [16] : λ u _r = 10 ^{- 12} и λ u _φ = 10 ^{- 13} .

Таблица 1 и рисунок 2 фиксируют диагностический признак некорректности дискретной задачи 1-го рода. В таблице через σ k обозначено k -е сингулярное число нормированной матрицы, σ 1 — максимальное сингулярное число, σ min — минимальное. Всего матрица имеет 18 сингулярных чисел, по числу столбцов. В УСР используется рабочий ранг r = 14 , выбираемый по спектральному разрыву и устойчивости решения. Численный SVD-ранг, определяемый по машинному порогу малости сингулярных чисел, равен 15. Различие обусловлено разными критериями определения этих рангов. Рабочий ранг выбирается на основе анализа спектрального разрыва и критерия устойчивости (см. подраздел 4.2), вследствие этого он, как правило, меньше; оставшиеся сингулярные числа подавляются при регуляризации. Величины j _min /j i = 1.51 • 10 ^{- 16} и cond( A eq ) = 6.61 • 10 ¹⁵ , где A eq — уравновешенная матрица системы, полученная нормировкой столбцов A на их евклидовы нормы (разброс исходных норм см. в Табл. 1) , а cond = σ ₁ /σ _min — её спектральное число обусловленности, следует рассматривать как численные индикаторы вырожденности нормированной матрицы в пределах машинной точности.

Таблица 1. Показатели дискретной системы для эллиптической области

Показатель	Значение Описание
Размер системы	2 Q х N =640х 18  640 вещественных уравнений коллокации и 18 неизвестных коэффициентов разложения
Число сингулярных чисел	N σ =18 равно числу столбцов, поскольку 2 Q > N
Рабочий ранг УСР	r =14 число сингулярных компонент, сохранённых после усечения
Численный SVD-ранг	15 число сингулярных чисел из 18, превышающих машинный порог
Отношение минимального и максимального сингулярных чисел	^ min / ^ i = 1.51 • 10 16 индикатор численной вырожденности матрицы на уровне машинной точности
Разброс норм столбцов	2.115 • 10 3   отношение максимальной и минимальной евклидовых норм столбцов исходной матрицы

Рис. 2. Сингулярные числа нормированной системы для эллиптической области; пунктирной линией показан порог рабочего усечения 10 ^-11 для отношения σ k /σ 1 ; точечной линией выделен рабочий ранг r = 14

• 5.2.    Прямоугольная область со смешанными граничными условиями

Во втором численном эксперименте рассматривалась плоская задача для ортотропного прямоугольника:

Г2 =[0,a]x [0,b],     a = 2 м , b =1 м .

Граница включала четыре участка:

дп=rL urR иГв и Гт, rL = {(0,x2):0 C X2 C b},

Г в = {(x i ,0):0 <  x i <  a},

Г я = {(a,x 2 ):0 C x 2 C b}, Г т = {(x _i ,b) :0 C x _i C a}.

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

^U | r L = 0 ,      t l r _B и Г т = 0 ,

t | _r R = ( t 0 ^{(2 x} 2 ^{/b —} 1) ^{, 0} ) ,

где t ₀ = 1 МПа. Таким образом, левый торец (индекс L ) был защемлён, правый торец R находился под нагрузкой, а две горизонтальные стороны (нижняя B и верхняя T ) были свободны от нагрузок.

Материалом служил углепластик T300/5208 [26] с безразмерными упругими постоянными, нормированными на модуль упругости E 2 для плоского напряжённого состояния. Матрица C в форме раздела 2 имела коэффициенты:

C ii = 17.65, C 22 = 1.00, C 12 = 0.28, С бб = 0.70.

В полной смешанной постановке неизвестными являлись напряжения на защемлённой стороне и перемещения на трёх остальных сторонах:

ti(0,X2), t2(0,X2), Ui(a,X2), U2(a,X2), Ui(xi,0), U2(xi,0), ui(xi,b), U2(xi,b), то есть восемь граничных функций. При расчётах использовалась симметрийная редукция относительно средней линии x2 = b/2:

u i (x i ,b-X 2 ) = -U i (x i ,X 2 ),     U 2 (x i ,b — x 2 ) = U 2 (x i ,X 2 ).

Отсюда ui(xi,b) = -Ui(xi,0),     U2(xi,b)= U2(xi,0).

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

t i (0,x 2 ), t 2 (0,x 2 ), u i (a,x 2 ), u 2 (a,x 2 ), u i (x i ,0), u 2 (x i ,0).

Параметры дискретизации прямоугольной области сведены в таблицу 2.

В углах типа «заделка–свободная граница» компоненты вектора напряжений имеют локальную степенную особенность, характерную для клиновидных областей в анизотропной теории упругости [3, 17, 18] . Анализ задачи для ортотропного клина дал значение показателя y = v — 1 = —0.15 .

Таблица 2. Показатели дискретной системы для прямоугольной области
Показатель	Значение	Описание
Независимые граничные функции	6	После симметрийной редукции верхняя свободная сторона восстанавливается по нижней
Число коэффициентов на одну функцию	N c =40	Используются полиномы Чебышёва 1-го рода T n по локальной координате, приведённой к [-1,1]; число членов выбирается по стабилизации граничных следов и рабочего ранга УСР
Коэффициенты гладкой части	6 N c =240	Число коэффициентов до наложения условий согласования в угловых точках
Независимые коэффициенты после согласования	N =234	Условия согласования на смежных сторонах уменьшают число свободных коэффициентов
Параметры коллокации	N ξ = N η =48	Число значений координатных параметров характеристического множества.
Комплексные коллокационные точки	Q =4 N ξ +4 N η =384	Для каждого фиксированного параметра учитываются четыре корня характеристического уравнения
Строки вещественной системы	2 Q =768	Получаются после разделения комплексной системы на действительную и мнимую части
Размер дискретной матрицы	2 Q X N =768x234	Число строк равно 2 Q , число столбцов — числу независимых коэффициентов
Число сингулярных чисел	N σ =min(2 Q , N )=234	Значение используется при анализе сингулярного спектра до регуляризации

Длина локального участка S t выбиралась равной S t /b = 0.006 . В регулярной части t ₂ (0,x ₂ ) использовалась регуляризация с оператором второй производной L = D ² по локальной координате, как и в подразделе 5.1; параметр A _t2 = 10 - 4 определялся автоматически по критерию GCV [16] . Остальные блоки стабилизировались усечением сингулярного спектра по критерию раздела 4. Рабочий ранг итоговой регуляризованной системы УСР составил r = 199 .

Результаты для сравнения с МКЭ при А _ш =0.30 приведены на рисунке 3. Контрольное МКЭ-решение получено на сетке из 60 х 120 четырёхузловых четырёхугольных элементов с билинейной аппроксимацией перемещений (в англоязычной литературе — Q4), что соответствует 14 762 степеням свободы. Достаточность сетки проверялась по стабилизации граничных L 2 -норм при переходе к более густым прямоугольным разбиениям. Локальный сингулярный элемент снижает ошибку в определении t ₂ (0,x ₂ ) ; расчёты остальных компонент демонстрируют малые относительные отклонения.

Расчёты показали, что при А _ш = 0.30 включение локальных сингулярных элементов снижает среднюю ошибку АГИУ относительно МКЭ с 3.261% до 1.021% , а ошибку t ₂ (0,x ₂) — с 16.929% до 2.108% (Табл. 3) . Таким образом, угловые точки не препятствуют применению дискретной АГИУ-схемы, но требуют согласованного с геометрией локального учёта особенностей поля напряжений.

Таблица 3. Относительные L 2 -ошибки АГИУ относительно МКЭ для прямоугольной области при А _ш =0 . 30

Компонента	Ошибка, %	Компонента	Ошибка, %
U 1 ( a,x z )	0.698	U 2 ( a,X 2 )	0.058
U 1 ( x 1 ,b )	0.184	U 2 ( X 1 ,b )	0.424
t 1 (0 ,X 2 )	2.656	t 2 (0 ,X 2 )	2.108

• 5.3.    Резонансы прямоугольной области

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

u l r L = 0 ,      u l r R = 0 ,      t l r B = 0 ,      t l r T = 0 .                              ⁽²⁶⁾

Однородность условий (26) означала отсутствие внешнего возбуждения. Собственная частота соответствовала такому значению λ ω , при котором задача (2) , (26) имеет нетривиальное решение. Эталонные частоты получались путём независимого решения прямой спектральной задачи на последовательности сгущающихся конечно-элементных сеток.

В силу симметрии области, материала и условий (26) относительно средних линий x 1 = a/2 и x ₂ = b/2 (при переносе начала координат в центр области — отражений x 1 н —x 1 , x ₂ н —x ₂ ) спектральная задача раскладывалась по четырём классам чётности (от A до D ), которые введены в подразделе 4.3. В этих подпространствах открылась возможность раздельного поиска собственных форм. Пробное пространство каждого расчёта при этом уменьшалось вчетверо, а близкие резонансы разных классов не сливались в один пик индикатора

Рис. 3. Граничные функции прямоугольной области при λ ω = 0 . 30 : АГИУ (штриховые линии с кружками) и МКЭ (сплошные линии)

и надёжно различались (Табл. 4, столбец «Класс»). Объединение результатов по четырём классам дало полный спектр рассматриваемой смешанной задачи, так что собственные частоты при таком разделении не терялись.

В МЧР (см. подраздел 4.2 [22, 23] ), где приближённое поле представляет собой комбинацию решений, точно удовлетворяющих уравнениям движения (2) , а собственная частота отвечает минимуму нормы граничной невязки, для условий (26) используется индикатор (24) . В данном расчёте бралось по 14 значений декартовых параметров N ξ = N η = 14 ; для каждого фиксированного параметра учитывались четыре комплексных корня. Пробное пространство содержало N = 112 комплексных функций, что после разделения на действительную и мнимую части соответствовало 224 вещественным столбцам. Невязка условий (26) аппроксимировалась по 100 точкам на каждой стороне; нормировка в знаменателе (24) вычислялась по сетке с количеством 28 х 16 внутренних точек. Начальная сетка по λ _ω выбиралась с учётом предварительной локализации резонансов, локальное уточнение проводилось в окрестности шириной 0.06 по λ ω . МКЭ-эталонное решение получалось на последовательности прямоугольных Q4-сеток от 80 х 40 до 240 х 120 элементов с последующей экстраполяцией Ричардсона по параметру сетки h — характерному размеру элемента. Погрешность билинейных элементов по частоте убывала как O(h ² ) , что позволяло по последовательности сгущающихся сеток оценить предельное значение А _ш при h ^ 0 . В таблице 4 через ε rel обозначено относительное отклонение результатов МЧР от МКЭ-эталона.

Таблица 4. Первые пять собственных частот прямоугольной области для условий (26) в расчётах МЧР и МКЭ

№	Класс	λ ω МЧР	λ ω МКЭ	ε rel , %
1	D	1.24127	1.24194	0.0538
2	C	2.47834	2.47982	0.0594
3	A	3.33697	3.33740	0.0129
4	D	3.78171	3.78403	0.0613
5	B	3.86388	3.86606	0.0565

5.4. Резонансы эллиптической области

Далее представлены результаты определения резонансных значений λ ω для области вида ортотропного эллипса с указанными выше константами и свободной границей:

t ¹ \дЯ = t 2 \ dn = ^{0 .}

Безразмерная частота задавалась параметром А _ш = wR 0 у/ p/E 2 , где, как и в подразделе 5.1, R 0 = ab = = 0.1 м.

В силу симметрии области и свободной границы спектральная задача редуцировалась к задаче на четверти эллипса. Собственные частоты определялись по индикатору (24) и сравнивались с МКЭ-решением обобщённой задачи на собственные значения. Исходное пробное пространство содержало 32 вещественные функции. Норма граничной невязки аппроксимировалась по 160 точкам на редуцированной четверти границы. Локализация каждого максимума Ф(А _Ш ) выполнялась в окрестности полушириной 0.08 по А _ш . Для каждого из четырёх резонансов начальный поиск проводился на равномерной сетке из 81 значения безразмерной частоты λ ω в интервале от (А МКЭ — 0.08) до (А МКЭ + 0.08) , после чего положение максимума индикатора уточнялось. МКЭ-эталон получался на сетке из 65 536 шестиузловых треугольных элементов с квадратичной аппроксимацией перемещений, число степеней свободы составляло 263 170.

В таблице 5 приведены данные вычислений на основе АГИУ и МКЭ. Через ε _rel обозначено относительное отклонение от МКЭ-эталона.

Таблица 5. Собственные частоты эллипса: АГИУ и МКЭ

№	Класс	λ ω АГИУ	λ ω МКЭ	ε rel , %
1	C	1.71868	1.71870	0.001254
2	D	1.87834	1.87836	0.001256
3	A	2.16455	2.16458	0.001255
4	B	2.42005	2.42008	0.001258

На рисунке 4 по оси ординат отложено значение индикатора Ф(А _Ш ) из (24) . Индикатор вычисляется отдельно для каждого из четырёх классов чётности (пробное пространство состоит из полей соответствующей чётности, поэтому каждая кривая отвечает своему классу). Найденный максимум сразу относится к определённому классу симметрии (см. Табл. 5, столбец «Класс»). Собственная частота в АГИУ определяется как абсцисса локального максимума кривой Ф(А _Ш ) своего класса.

Рис. 4. Индикатор Ф ( Х ^ ) для первых четырёх собственных частот эллиптической области, рассчитанный на основе АГИУ (сплошные линии), и значения λ ω , вычисленные МКЭ (пунктирные линии)

• 6.    Заключение

Построена и численно реализована схема альтернативных граничных интегральных уравнений (АГИУ) 1-го рода с гладкими ядрами для установившихся колебаний ортотропного тела в плоской постановке. Система АГИУ формируется через условия разрешимости на характеристическом множестве дифференциального оператора и не требует явного построения фундаментальных и сингулярных решений. Предложенная схема приводит к дискретизации только граничных неизвестных и применима как к прямым смешанным постановкам, так и к спектральной задаче о собственных частотах.

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