CC BY
6
УДК 519.635.1+519.632.4+539.3 Б01: 10.14529/шшр230303
РАЗРАБОТКА И ВЕРИФИКАЦИЯ УПРОЩЕННОГО НР-ВАРИАНТА МЕТОДА КОЛЛОКАЦИИ И НАИМЕНЬШИХ КВАДРАТОВ ДЛЯ НЕРЕГУЛЯРНЫХ ОБЛАСТЕЙ
Л.С. Брындин1'2, В.А. Беляев1, В.П. Шапеев1
Институт теоретической и прикладной механики им. С.А. Христиановича СО РАН, г. Новосибирск, Российская Федерация
2Новосибирский государственный университет, г. Новосибирск, Российская Федерация
Предложен, реализован и верифицирован новый высокоточный Ир-вариант метода коллокации и наименьших квадратов (Ър-МКНК) численного решения эллиптических задач в нерегулярных областях. При построении приближенного решения использовались граничные нерегулярные ячейки (н-ячейки), отсеченные границей области от ячеек прямоугольной сетки, и их законтурные части для записи уравнений коллокации и условий согласования. В малых и (или) вытянутых несамостоятельных н-ячейках отдельное решение не строилось, а продолжалось из соседних самостоятельных ячеек, в которых использовалась внешняя (и внутренняя в многосвязной области) часть границы области, заключенная в этих несамостоятельных н-ячейках, для записи краевых условий. Такой подход существенно упростил компьютерную реализацию разработанного Ир-МКНК по сравнению с предыдущим хорошо зарекомендовавшим его вариантом, не потеряв при этом своей эффективности. Показана возможность уменьшения степени переопределения системы линейных алгебраических уравнений по сравнению с ее значениями в традиционных вариантах МКНК при решении бигармонического уравнения. Проведено сравнение с результатами других работ с демонстрацией преимуществ нового подхода. Приведены результаты расчетов кольцевых пластин различной толщины на изгиб в рамках теорий Кирхгофа - Лява и Рейсснера - Миндлина с помощью Ир-МКНК, демонстрирующего отсутствие сдвигового запирания.
Ключевые слова: метод коллокации и наименьших квадратов; теория Кирхгофа-Лява; теория Рейсснера-Миндлина; бигармоническое уравнение; нерегулярная область.
Введение и постановка краевых задач
Эллиптические задачи встречаются во многих разделах механики сплошных сред. Примером одной из них является расчет изгиба пластин - распространенных несущих элементов многих конструкций. Анализ их напряженно-деформированного состояния (НДС) при изгибе вызывает ряд сложностей из-за различных краевых условий закреплений и геометрий расчетных областей [1]. Указанные обстоятельства наряду с рассмотрением различных теорий описания НДС пластин [2-4] усложняют строгий математический анализ и поиск точных решений задач, зачастую являющихся еще плохо обусловленными и с малыми параметрами при старших производных. В связи с этим для нахождения их решений необходимо разрабатывать и верифицировать численные методы, способные решать подобные проблемы. В данной статье для задачи статического изгиба однослойных изотропных пластин рассмотрены теории Кирхгофа - Лява (ТКЛ) [2,3] и Рейсснера - Миндлина (ТРМ) [1,3], в рамках которых требуется рассчитать их НДС с помощью эффективного численного метода.
Напомним, что расчет изотропных пластин с внешней границей y0 и внутренними границами Yi (i = 1, • • •, N, NY - количество отверствий) на изгиб в рамках ТКЛ сводится к решению бигармонического уравнения в области П С R2 [2]
d4w д4w д4w q . ,
__L 9__i__= — Г Л
дх\ дх\дх22 дх4 D'
где w(x 1,^2) - прогиб, q(x 1,^2) ~~ поперечная нагрузка, D = ———■—— - жесткость
12(1 — V )
пластины при изгибе, h - толщина пластины, E - модуль Юнга, v - коэффициент Пуассона.
Для уравнения (1) рассмотрим задачу с краевымии условиями на границе области, которые в случае защемленного края имеют вид
dw , ,
w = 0, — = 0, 2
дп
а в случае шарнирно закрепленного края
w = 0, Mn = MX1 cos2 a + Mx2 sin2 a + 2MX1X2 sin a cos a = 0, (3)
где n = (cos(a), sin(a)) - внешняя единичная нормаль к соответствующей границе области, a - угол между п и осью xi, отсчитываемый от x\ против хода часовой
тл (д2w д2w\ (д2w d2w\
стрелки, моменты МХ1 = -ü I + ищ \ , МХ2 = -ü I + J > МХ1Х2 =
д2w
^^ 8x18x2'
В последние два десятилетия исследователями ведутся активные исследования по разработке и применению метода конечных разностей (МКР) [5,6], метода конечных элементов (МКЭ) [7], спектральных методов [8] и других для решения краевых задач для бигармонического уравнения (1) в нерегулярных областях. Подробный обзор по современному состоянию численных методов их решения приведен в нашей свежей статье [9] и цитируемых там работах. В ней же продемонстрированы преимущества и достоинства реализованных h-, p- и hp-вариантов метода коллокации и наименьших квадратов (h-, p- и hp-МКНК) решения (1) в областях со сложной геометрией с краевыми условиями, содержащими производные низкого и высокого порядков, в сравнении с результатами, приведенными в [5,7,8] и других работах.
Однако компьютерная реализация предложенных в [9] подходов для построения решения в нерегулярных областях относительно сложна. В сеточных вариантах МКНК исходная нерегулярная область включалась в прямоугольник, который покрывался регулярной сеткой с прямоугольными ячейками. Малые и (или) вытянутые несамостоятельные нерегулярные ячейки (н-ячейки) присоединялись к соседним самостоятельным ячейкам по следующему правилу: несамостоятельная н-ячейка присоединялась к той соседней самостоятельной ячейке, с которой она имела наибольшую длину стороны и (или) части стороны, расположенной внутри области, среди всех ее сторон, общих с другими соседними самостоятельными ячейками. При этом программирование такого подхода требовало больших трудозатрат, так как возникала необходимость вводить много различных типов ячеек с учетом их ориентации в двумерном пространстве. Кроме того, в случае многосвязной области в [9] изменение формы отверстий приводило бы к заметным изменениям в программном коде.
Расчет изгиба пластин в рамках ТРМ, учитывающей поперечные касательные напряжения, сводится к решению системы дифференциальных уравнений, состоящей
из трех уравнений с частными производными (УЧП) второго порядка [3]. В случае изгиба однослойной изотропной пластины разрешающая система для определения прогиба и>(ж1,ж2) и углов поворота нормали (х, х2) и (х, х2) срединной поверхности вокруг осей х2 и х1 соответственно имеет следующий вид:
-КА^ - К А^ - К Л
дх2
дх2
дх1
- К А
д2
т
дх1
(£>12 + Азб)п - КА(рХ1 - К А— + + В66-
(^12 + Д66)
дх1 дх2 ГХ1 д^ ' 11 дх1
^ -КА^-КА^ + Вп9^
дх2
Х2
дх1 дх2
дх2
дх2
+ В
д2^
66"
Х2
дх1
0,
0,
где КС = 5/6 - сдвиговой коэффициент Тимошенко, А = СЛ, В12 = Ек*^, В66
Л3 Е
(4)
С/?3 12 '
Яп = Ек*, к*
С
12(1 - V2)' ~ 2(1 + V)'
В этом случае по сравнению с ТКЛ в (4) при старших производных и возникают малые параметры (см. раздел 1.3 [10]). Также существенным затруднением численного решения этой задачи МКЭ является наличие эффекта сдвигового запирания: по мере уменьшения толщин приближенное решение все хуже удовлетворяет требованиям, допускающим нулевые деформации поперечного сдвига [11, 12]. Существуют разные способы преодоления этого обстоятельства: смешанная аппроксимация [12], повышение степени полиномиальных приближений [11,12] и др. Так, в [11] отмечается, что одним из эффективных подходов в МКЭ является использование конечных элементов 3-й степени или выше (в зависимости от толщины пластины). Однако, к сожалению, далеко не всегда удается полностью избежать этого эффекта, а получается только уменьшить его влияние [12].
Впервые спектральный МКНК для расчета изгиба прямоугольной пластины в рамках ТРМ был реализован и применен в [4]. При этом наличие существенных особенностей (разрывы производных, точки ветвления) не позволяет эффективно применить «в лоб> многие спектральные численные методы [4,8-10], в которых приближенное решение строится в виде единого полинома достаточно высокой степени во всем прямоугольнике, включающем исходную нерегулярную область, даже если особенность находится внутри отверстия [9]. Одним из способов решения этой проблемы в многосвязных областях является построение и улучшение сеточных вариантов численных методов, чему и посвящена данная работа.
В [1] в рамках ТРМ даны аналитические решения задачи изгиба шарнирно-закрепленной круглой пластины с условиями [3]:
т
0, Мп = 0, ^ = сов(а) - йт(а) = 0,
(5)
где МЛ1
В
д^;
11"
Х1
В
66
+
д^;
дх1
+ В
д^
12'
Х2
дх2
Мх2 =
В
д^
12'
Х1
дх1
+В
д^;
11"
Х2
дх2
Мх
Х2
дх2 дх1
Точные решения являются важными для верификации численных методов в нерегулярных областях. Так, в [13] численное решение изгиба защемленной круглой пластины в рамках ТРМ сравнивалось с точным решением из [1]. К сожалению, авторы [13] не привели сведений об исследовании точности и сходимости численных методов при различной дискретизации, ограничившись рисунками приближенного и точного ее решений.
Настоящая работа направлена на существенное упрощение Ьр-МКНК с сохранением его эффективности, верификацию и сравнение с результатами других работ. Также отметим, что здесь он впервые применен для расчета НДС пластин непрямоугольных форм в рамках ТРМ.
1. Ьр-МКНК
Опишем новый Ьр-МКНК на примере решения (1), (2). Область П поместим в прямоугольник (рис. 1 а) и покроем его прямоугольной сеткой размера N1 х (рис. 1 б). Напомним, что граничными н-ячейками считаются ячейки, имеющие непустое пересечение с границами области и не совпадающие с исходными материнскими прямоугольными ячейкам, которые содержат их внутри себя (см. рис. 1 б, ячейки 1-8) [9]. Если центр материнской ячейки находится внутри П, то такие ячейки являются самостоятельными (см. рис. 1 б, ячейки 1-3, 8), иначе - несамостоятельными (см. рис. 1 б, ячейки 4-7). В каждой ]-й самостоятельной ячейке, ] = 1,..., Мсец3, где Мсец3 - их количество, введем локальную систему координат
У1
X1 X ij
У 2
X2 — X2 j
ho.
(6)
где (х^,х2]) - центр ячейки, 2Н1, 2Н2 - ее размеры вдоль осей Ох1 и Ох2 соответственно.
о Sr ° ' / 8 [ О О ! : о о 1 Q : о о : : о 3 о > : о о и'З
Ф / 6 О 0 ! 9 О О ! : о 2 о : о о \7/ >- / о /то
<j> О \ О О О 4/ о
аб
Рис. 1. Область из примера 1 (а) и фрагмент сетки (б) при К = 2 и М1 = М2 = 10, где значок о обозначает точки коллокации, х - точки согласования, □ - точки записи краевых условий, о - центры несамостоятельных н-ячеек
Приближенное решение в каждой ]-й самостоятельной ячейке ищем в виде [9]
к к-i1
Whj (У1,У2) = Y1 Ciii ¿1=0 i2=0
ili2,Jy11 y22 ,
(7)
где вг1 г2,] - неизвестные коэффициенты. Для их определения в каждой самостоятельной ячейке выпишем переопределенную локальную СЛАУ из уравнений коллокации, согласования и краевых условий (для граничных ячеек).
Уравнения коллокации - результат аппроксимации искомым решением (7) уравнения (1), записанного в локальных переменных (6) [9]:
, 1 д4 и>^__2 д4 и>^ 1 д4 и>^\ _
Рс[ ду\ + тду1ду1 Ъ ду4 ^
где рс > 0 - весовой множитель уравнения коллокации. Их выпишем в точках кол-локации, равномерно распределенных в самостоятельных ячейках (рис. 1 б) в количестве Мс = (К — Мр)2, Мр - параметр, регулирующий степень переопределения СЛАУ. При этом небольшая часть точек коллокации в граничных ячейках могла располагаться вне области (см. рис. 1 б, ячейка 8). В этом случае предполагается, что уравнение справедливо также и в малой окрестности ее границы.
Условия согласования обеспечивают непрерывность решения и его производных
[9]:
дцуц _ . ( дтн
Рто ^М I Ртх п Рто ^И \ Рт\ п )
дп7- дп7
7 7 (9)
Рп12 о ~г Рпг 3 ^ -3 Рп12 о ~г Рпг з ^ -з ) дп:- дп3 оп- дп3
где йь - приближенное решение в соседней ячейке, п7 - внешняя нормаль к границе ]-й ячейки, рто ,рт1 ,рт2 ,ртз > 0. Условия (9) записывались в Мт = (К — Мр) точках, равномерно распределенных на общей стороне двух самостоятельных ячеек. При этом в некоторых случаях они могли быть также выписаны вне области (рис. 1 б, ячейка
3).
Замечание 1. Покажем, что выполнение первого уравнения (9) влечет непрерывность отдельно функции и ее первой производной по нормали. Для этого рассмотрим условие согласования между соседними ячейками с номерами ] и ] + 1. Выпишем это условие для каждой из них с учетом того, что п7+1 = — п7 (следовательно,
—- = — ——), так как внешние нормали к оощеи стороне у двух соседних ячеек
дп]+1 дп]
противоположно направлены:
Рто^з + Рпг - = + Ргщ—^-, (Ю)
дйнз+1 дй^^
Рт0 «'/гз+1~Рт1 - = Рт0 «'/гз ~ Рт1 Т,-• (11)
дп? д п.,
Попарно складывая и вычитая уравнения (10) и (11), получим
дгуН] _ дгуы+х дп7- дп.,-
2рто йщ = 2рто йЬ]+1-
дйь] дйь]+1 . „
Следовательно, и'м = 1 и —— = —--в любой рассматриваемой точке согла-
дп] дп]
сования. Аналогичные рассуждения справедливы для второго уравнения (9) с учетом
д2 д2 д3 д3
т—^з— = —--т и ——г— = — ——т. Отметим, что использование симметричных условий дп2+1 дп2 дп3+1 дп3
согласования (9) уменьшает общее количество уравнений в переопределенной СЛАУ
и, как правило, ускоряет сходимость итерационного процесса по сравнению с вариантом раздельной записи условий непрерывности функций и независимо от них -условий непрерывности производных [14].
Краевые условия, домноженные на весовые множители рь0,рь1 > 0, имели вид [9]
Pbi
Pbo Wh dwh
dn
0, 0.
(12)
Для расстановки Nb точек записи краевых условий определялась часть границы области (как в случае внешней y0, так и внутренних Yi, i = 1,... , NY, в случае много-связности), заключенная в прямоугольник, найденный по следующему правилу. Если у j-й самостоятельной ячейки примыкающая к ней слева ячейка является несамостоятельной, прямоугольник выбирается объединением материнской прямоугольной j-й ячейки и половиной соседней, отсеченной вертикальной линией, проходящей через центр этой несамостоятельной н-ячейки. Так, на рис. 1 б для записи краевых условий в ячейке 9 выбрана часть границы области, заключенная между двумя вертикальными линиями, одна из которых пунктирная и проведена через центр ячейки 6, а другая является правой стороной ячейки 9, и горизонтальными линиями, ограничивающими ячейку 9. Если левый сосед отсутствует для j-й ячейки, то в качестве прямоугольника выбирается материнская j -я ячейка.
Аналогично рассматриваются также и ситуации, когда необходимо учесть границы Yi, i = 0,... , Ny, для записи краевых условий в самостоятельной ячейке, к которой примыкают несамостоятельные н-ячейки справа, снизу и сверху, и которых может быть одновременно несколько. Например, на рис. 1 б для записи краевых условий в ячейке 1 выбрана часть границы области, заключенная между двумя пунктирными линиями, проходящими через центры ячеек 4 и 5, и горизонтальными линиями, ограничивающими ячейку 1.
В итоге на найденной части границы точно и равномерно располагались Nb точек. Величина Nb равнялась произведению (K — Np) на количество отсутствующих и несамостоятельных н-ячеек, примыкающих к текущей самостоятельной снизу, справа, сверху и слева.
Решение в каждой точке несамостоятельной н-ячейки продолжим из самостоятельной ячейки, расстояние от центра которой до рассматриваемой точки наименьшее по сравнению с другими самостоятельными ячейками. Например, для определения решения в ячейке 7 в ее темно-серой части решение продолжается из ячейки 2, а в светло-серой - из ячейки 3, как показано стрелками на рис. 1 б.
Замечание 2. Значения весовых множителей pc,pm0 ,pmi ,pm2 ,pm3,Pb0,Pb1 в МКНК при применении взвешенного метода наименьших квадратов влияют на обусловленность локальных СЛАУ, скорость сходимости итерационного процесса и точность решения задачи. Для каждой конкретной задачи существуют целые области их значений, которые обеспечивают сходимость их решений и находятся в численных экспериментах (см., например, [15, 16] и цитируемую там литературу). Здесь мы для задания весов воспользовались опытом предыдущих исследований [9,15].
Глобальная СЛАУ, являющаяся объединением всех локальных, решалась с помощью метода итераций по подобластям [17]. Итерационный процесс прекращался, когда
max K+L-— cSii2,j1 (13)
il ,i2,J 1 2'J 1
где ^, с|+21 ^ - коэффициенты в представлении (7) на 5-й и (в + 1)-й итерациях, в = 0,1,... , Мцег, Мцег - количество итераций, е - малая величина.
Замечание 3. При решении системы (4), состоящей из УЧП второго порядка, каждая неизвестная функция аппроксимировалась выражениями вида (7). В качестве условий согласования выписывалось только первое уравнение (9) для каждой неизвестной. Величины N5 определялись по указанным выше формулам с ис-
пользованием в них в качестве К самой старшей степени из всех рассматриваемых полиномов.
Для улучшения сходимости и существенного сокращения количества итераций и времени расчетов в методе итераций по подобластям применялись: метод ускорения итерационного процесса, основанный на подпространствах Крылова [16]; операция продолжения для задания начального приближения на подробной сетке при переходе на нее с более грубой [9]; распараллеливание программы, написанной на языке С+—Н, с помощью ОрепМР в комбинации с обходом ячеек сетки с применением красно-черного упорядочивания; свойство специального вида матриц и правых частей локальных СЛАУ в МКНК [17].
2. Результаты численных экспериментов и их обсуждение
В представленных ниже численных результатах абсолютная и относительная погрешности определялись соответственно по формулам
= max ( max |w(Xim,X2m ) - Whj (xim ,X2m)|),
j= l,...,n1n2 m=1,...,100
EW
max ( max |w(Xim,X2m) - Whj (xim ,X2m )|)
j=l,...,Nin2 m=l,...,100
max ( max |w(xim, x2m)|)
(14)
где (х1т, х2т) - равномерно распределенные точки в каждой материнской ячейке для подсчета в них погрешности. Если (х1т, х2т) / П, то в ней погрешность не рассматривалась. Величины (14) аналогично определялись для , и напряжений, которые в ТКЛ имеют вид [2]
(71 = -e(wxlxl + uwxnxn ) 1
Ехз
E
1 - V
2
в случае ТРМ [3]
(Ti = Е i
xi
+ Е
Х2
12'
dx2
. dw
02
E
-e(wX2X2 + uwXlXl Ехз
1 + v
), ^12 Хз G [-h/2, h/2],
(72 = E
xi
12'
X2
dx2
. dw
0"23 = Cr ( — + (/?.T2
E,
ai2 = Gx3 Ехз
1 - v
2
—Еш
XiX2 )
<9</?.Tl <9</?.T2
dx2
E
12
dx, Evx3 1 -//2
Порядок сходимости погрешности приближенного решения определялся как R =
l°g2
E
N/2
E
N
EN/2 и EN - погрешности на сетках размера N,/2xN2/2 и N,xN2 соот-
ветственно. Расчеты были проведены на персональном компьютере 12th Gen Intel(R) Core(TM) Í7-12700H 2,70 GHz, DIMM DDR5 2400 MHz 16 Gb. Для распараллеливания вычислительной программы с помощью OpenMP использовалось 20 потоков.
Далее через г801 обозначено время выполнения итерационного процесса в секундах. В расчетах на самой грубой сетке полагалось с0^ ^ = 0, 4.
Пример 1. В этом примере, косвенно относящемуся к механике деформируемого твердого тела, одной из целей являлось сравнение нового Ьр-МКНК с предыдущим вариантом МКНК [9] и МКР [5,6] на решении тестовых краевых задач Дирихле для бигармонического уравнения. Для этого заранее известное точное решение задачи подставлялось в (1), (2) и вычислялись их правые части.
Рассмотрим сначала предложенную в работе [5] задачу (1), (2) с тестовым решением , х2) = х12 + х22 + еХ1 сов(х2) в области в форме звезды (см. рис. 1 а), внешняя граница которой задана следующими соотношениями:
х1 = (0, 6 + 0, 25 8т(5£)) соэ(^), г е [0, 2п), х2 = (0, 6 + 0, 25 яп^)) в1п(г).
В табл. 1 приведены результаты расчетов, полученных с помощью разработанного здесь Ьр-МКНК. В ней Мед обозначает количество уравнений в переопределенной локальной СЛАУ. В [9,17] величина Мед = 48 при К = 4 и Мед = 84 при К = 6. Здесь, как ив [9], при применении взвешенного метода наименьших квадратов были взяты следующие веса: рс = Л? , рто = рЪо = Рь1 = 1, Рт 1 =0, 01, рт2 = Л^, Рт3 = 0, 015Л1Л2. Для остановки итерационного процесса в (13) при К = 4 использовалась величина е = 10-10 на всех сетках, кроме 80x80, где она равнялась 10-11. При К = 6 величина е = 10-12 на сетках 5x5, 10x10 и 20x20, а на сетке 40x40 она была взята равной 0, 5 • 10-13. Размер базиса подпространства Крылова в применяемом здесь методе ускорения итерационного процесса был равен 40 [16].
Таблица 1
Результаты численных экспериментов в примере 1
К = 4 ( Neq = 33, Np = 1) К = 4 ( Neq = 48, Np = 0)
NiXN2 К R Niter t sol К R Niter t'Sol
10x10 8,60E-6 — 41 0,015 8,40E-6 — 41 0,016
20x20 l,31E-6 2,71 81 0,047 l,09E-6 2,94 81 0,062
40x40 7Д5Е-8 4,19 81 0,172 6,41E-8 4,08 81 0,219
80x80 3,92E-9 4,18 107 0,750 3,84E-9 4,06 108 0,969
К = 6 ( Neq = 48, Np = 2) К = 6 ( Neq = 84, Np = 0)
NixN2 К R Niter t sol К R Niter t'Sol
5x5 ЗД5Е-7 — 81 0,015 9,20E-8 — 48 0,016
10x10 5,06E-9 5,95 71 0,031 4,38E-9 4,41 41 0,047
20x20 1,69E-10 4,90 81 0,125 1,44E-10 4,92 81 0,187
40x40 9,99E-12 4,08 81 0,453 5,42E-12 4,73 125 0,906
Результаты табл. 1 практически не отличаются по значениям Еа, Я и при одинаковых К и демонстрируют возможность существенного сокращения количества уравнений в переопределенных СЛАУ. Это позволяет в первом случае уменьшить время решения задачи в сравнении со вторым случаем. Максимальные числа обусловленности локальных СЛАУ ^тах в спектральной норме на рассмотренных здесь сетках изменялись от С1 • 102 до С2 • 103 и от С3 • 103 до С4 • 104 (С1, С2, С3, С4 - константы) при К = 4 и К = 6 соответственно как для случая, когда было уменьшение количества уравнений, так и для традиционного случая, когда сокращения не происходило. Причем различия ^тах при одинаковых К и N1 x N2 были незначительными.
Из табл. 1 также видно, что сочетание различных способов ускорения итерационного процесса решения СЛАУ позволяет на персональном компьютере решать УЧП четвертого порядка не только с высокой точностью, а также достаточно быстро и экономично.
В [5] при решении этой задачи МКР второго порядка аппроксимации Еа изменялась от 8,15Е-4 до 1,83Е-5 на сетках 64x64,..., 512x512. В [9] МКНК при К = 4 (7) удалось достигнуть точность 1,46Е-5, 1,70Е-6, 2,17Е-7 и 2,96Е-8 на последовательности измельчающихся вдвое сеток 10 х 10,..., 80x80. Таким образом, с помощью нового Ьр-МКНК удается достигнуть лучшие по точности результаты в сравнении с результатами [5] и [9], а в последнем случае еще и с помощью более простого алгоритма.
Реализация Ьр-МКНК для произвольной степени полиномов (7) автоматизирована, а также является достаточно простой для нерегулярных областей. Для других численных методов одновременный учет этих факторов может привести к серьезным трудностям. Кроме того, нам известно всего несколько работ по решению бигармо-нического уравнения в нерегулярных областях с помощью МКР (см. [5,6] и другие упомянутые работы во введении к статье [9]). Это наблюдение также подчеркивают и другие специалисты [8], относя сюда МКЭ помимо МКР. Отметим также недавнюю статью [6], в которой предложен МКР четвертого порядка аппроксимации для решения бигармонического уравнения в нерегулярных областях. Работ по построению МКР более высокого порядка для решения данной задачи мы в настоящий момент не встречали.
На рис. 2 приведены графики сходимости Ьр-МКНК и МКР [6], полученные при решении задачи Дирихле (1), (2) для бигармонического уравнения с тестовыми решениями и) = ех1+х2 (рис. 2 а) и = 10 ( х? вт^та^) Н--—) (рис. 2 б) в области из
V 1 + х\)
двух пересекающихся кругов, граница которой задана следующими соотношениями:
'хл = (0,4|сов(*)| + \/0, 25 - 0,16 8П12(£)) сов(£), I € [0, 2тг), = (0, 4| сов(*)| + а/0, 25 — 0,16вт2(г)) вт(г).
-4
-6
м
-10
-12
= 2,38
Д = З/28^Е
= 3,99 ^ \
10x10 20x20 40x40 NxN
а
80x80
10x10 20x20 40x40 NxN
б
80x80
Рис. 2. Графики сходимости численных решений на последовательности сеток в логарифмическом масштабе, где значок □ обозначает результаты, полученные МКНК при К = 4, △ - при К = 6, о - при К = 8, о - результаты, полученные МКР [6]
Сплошными линиями на рис. 2 показаны полиномиальные аппроксимации первой степени дискретно заданных значений порядков сходимости Я, полученных с помо-
щью метода наименьших квадратов, и указаны их осредненные значения. При этом для МКР они были взяты из работы [6]. При построении кривой на рис. 2 а не учитывались результаты на сетке 80x80 при K = 6 в МКНК, так как уже не происходило их улучшения. Также поступали и авторы [6], не учитывая результаты на сетке 60x60 на рис. 2 б.
Из представленных результатов видно, что в МКР в первом случае (см. рис. 2 а) понижается R по сравнению со вторым (см. рис. 2 б), чего не происходит в hp-МКНК. Повышение степени полиномов K приводит к заметному улучшению точности решений и порядку сходимости, позволяя получить решение задачи с более высокой точностью, чем МКР [6].
Пример 2. Рассмотрим кольцевые изотропные пластины (с параметрами E = 200 ГПа, v = 0, 28) различной толщины h с центральным отверстием xf + x2 = 12. Проведем расчет их НДС в рамках ТКЛ и ТРМ, когда они шарнирно закреплены по обоим краям и находятся под действием постоянной нагрузки q. В рамках ТКЛ точное решение краевой задачи (1), (3) имеет вид [2]
w = + Cir2(log г - 1) + Щ- + С3 log г + С4, г = Jxf + xl 64D 4 V
а в рамках ТРМ при решении (4), (5)
/logr r2logr r2 \ C6r2 qr2 qr4
(r Сч с, qr2 \ / cv C5 qr2 \
где константы Ci-C4 и C5-C8 определяются из решения СЛАУ размера 4x4 при учете краевых условий, соответствующих рассмотренным теориям пластин. В случае ТРМ для получения точного решения нами была применена техника из [1], где при
интегрировании уравнений слагаемые, содержащие logr и —, не отбрасывались.
r
В табл. 2 приведены результаты численных экспериментов, где S обозначает отношение диаметра внешней окружности xf + x2 = 52 к толщине пластины. Для ТРМ в скобках приведены степени аппроксимирующих полиномов K (7) последовательно для w, <рХ1 и ^Х2. Для трех случаев S = 10, 20,100 рассматривалась нагрузка q = 100, 10, 0, 01 МПа соответственно. Погрешность найденных напряжений рассматривалась при x3 = h/2.
Численные результаты, полученные hp-МКНК при расчете изгиба в ТКЛ при одинаковых NfxN2 и K для рассматриваемых S и q, практически не отличались по значениям погрешностей, поскольку в (8) изменялась только правая часть. Поэтому в табл. 2 в качестве примера приведены результаты для S = 20. Отметим также, что полученные численные результаты в этом примере практически не отличались, как и в примере 1, для Np = 0 и Np = 1, 2 при K = 4, 6 соответственно. Однако из-за ограничения на объем статьи они здесь не представлены.
В расчетах изгиба пластин в рамках ТРМ установлено, что при S = 10 и S = 20 отличия значений погрешностей были меньше по сравнению со случаем S = 100. Например, при использовании полиномов (4, 3, 3) при S = 100 на сетке 80x80 была достигнута точность E'W = 3,96E-4, а использование полиномов (6, 5, 5) при S = 20 на сетке 80x80 позволило достичь точность EW = 3,11E-7.
Из табл. 2 также видна высокая точность и порядок сходимости характеристик НДС, значения которых повышаются с увеличением степеней аппроксимирующих
Таблица 2
Результаты численных экспериментов в примере 2
ТКЛ, к = 4, мр = 1, 5 = 20 ТКЛ, к = 6, мр = 2, 5 = 20,
20x20 40x40 80x80 160x160 20x20 40x40 80x80 160x160
е? я 5,32Е-2 2,70 6,07Е-3 3,13 7,82Е-4 2,96 9,96Е-5 2,97 3,95Е-4 5,91 4,71Е-5 3,06 7,26Е-6 2,69 3,44Е-7 4,41
е.^ я 1,44Е-1 1,96 3,01Е-2 2,25 8,36Е-3 1,84 1Д8Е-3 2,82 1,70Е-2 3,32 1,32Е-3 3,68 1,97Е-4 2,74 1,00Е-5 4,30
я 1,54Е-3 1,85 3,97Е-2 1,95 8,09Е-3 2,29 1,60Е-3 2,33 2,30Е-2 3,12 2,54Е-3 3,17 2Д8Е-4 3,56 1,66Е-5 3,71
ТРМ, (4, 3, 3), АГР = 1, 5 = 10 ТРМ, (6, 5, 5), Жр = 1, 5 = 100
мгхм2 10x10 20x20 40x40 80x80 10x10 20x20 40x40 80x80
е? я 2,98Е-2 3,15Е-3 3,24 2,30Е-4 3,77 4,86Е-5 2,24 2,23Е-2 4,80Е-4 5,53 3,87Е-5 3,63 2,49Е-6 3,95
ег1 я 3,94Е-2 3,80Е-3 3,37 1,23Е-3 1,62 1,88Е-4 2,70 3,48Е-2 8,07Е-4 5,43 7,78Е-5 3,37 5Д7Е-6 3,91
е,^ я 1,67Е-1 2,60Е-2 2,68 1,03Е-2 1,33 1,57Е-3 2,71 8,70Е-2 2,97Е-3 4,87 3,30Е-4 3,16 2,78Е-5 3,56
е^ я 1,90Е-1 4,34Е-2 2,13 9,48Е-3 2,19 2Д8Е-3 2,12 6,62Е-2 8,91Е-3 2,89 8,73Е-4 3,35 7,95Е-5 3,45
е^ я 3,43Е-1 6,89Е-2 2,31 1,40Е-2 2,29 7,30Е-3 0,93 1,88Е+0 7,48Е-2 4,65 1,67Е-2 2,16 9,61Е-4 4,11
полиномов. Кроме того, представленные численные результаты подтверждают отсутствие эффекта сдвигового запирания в Ьр-МКНК для решения задачи изгиба тонкой пластины в рамках ТРМ.
Хорошо известно, что отличия между характеристиками НДС в ТКЛ и ТРМ будут увеличиваться при уменьшении Б [4,10]. На рис. 3 а приведены профили прогиба т пластины с толщиной к = 1 м (Б = 10) при х2 = 0, Х\ > 0, полученные МКНК в расчетах по ТКЛ и ТРМ. Максимальное относительное отклонение прогибов составляет приблизительно 20 %. Для круглой пластины аналогичной толщины относительная разница будет составлять 4 %.
На рис. 3 б приведено распределение напряжения а13 по толщине при Х1 = 1, х2 = 0. Здесь также для восстановления касательных напряжений применена методика из [3, раздел 5.2.2], по которой можно получить а13, проинтегрировав первое из уравнений равновесия трехмерной теории упругости по переменной х3:
хз
[ (да 1 0(712 \ ,
<713 = / ~--Ь — ах3.
} \дх1 0x2 /
2
Однако даже восстановленные а13 в ТКЛ отличаются более чем на 50 % от их значений, полученных в ТРМ при моделировании кольцевой пластины (см. рис. 3 б). В случае круглой пластины относительные отклонения для восстановленных а13 составляют уже 18 %.
Точным решением задачи изгиба шарнирно закрепленной круглой пластины без отверстия в ТКЛ (1), (3) [2] и ТРМ (4), (5) [1] для прогиба т в обоих случаях является полином четвертой степени, а в последнем случае для <рХ1 и (рХ2 — полином третьей
/ / ..... s \ .. \
/ /. / / / "4 \ \ Л
/7 /'/ \\ \\
// e' * r •A \ \
0.4 0.2 0.0 -0.2 -0.4
хз
! 4
N 4 4 \
J
!
—
1
3
х1
0
600
200 400
0-13 (1,0, х3)
а б
Рис. 3. Графики профиля прогиба ,ш(х1, 0) в расчетах по ТКЛ (• • •) и ТРМ (- • -) (а)
и напряжения а13(1, 0,х3) для ТКЛ (исходные (---), восстановленные (• • •)) и ТРМ
(исходные (--), восстановленные (- • -)) (б)
степени. Выбирая в МКНК в (7) степень К = 4 для прогиба и К = 3 для углов поворота, получаем решение с точностью порядка 10-14-10-15 на сетке 2x2.
Замечание 4. Порядок сходимости МКНК и правильность полученных результатов здесь и в других работах обосновываются только результатами численных экспериментов. Отметим, что аналогичный подход в МКНК был применен при решении линейных УЧП второго [10,14] и четвертого порядка [9,15,17] эллиптического типа, а также нелинейных УЧП [16,18], в том числе при исследовании Ьр-подхода с различными степенями полиномов К [9,10,17,18]. Можно считать, что апостериорный анализ в этом случае в некоторой степени выступает альтернативой строгим априорным оценкам, получение которых в настоящее время представляет собой весьма трудную, но актуальную задачу.
Заключение
Заметное упрощение реализации Ьр-МКНК, возможность с помощью него достигать повышенную точность решения задач в нерегулярных областях, его стереотипность как с точки зрения применения в случае различных постановок задач, так и использование при их решении произвольных степеней аппроксимирующих полиномов делают этот метод конкурентоспособным с популярными на сегодняшний день другими численными методами, прежде всего с МКЭ. При этом одновременный учет вышеуказанных факторов может быть достаточно трудоемким для разработки и реализации вариантов многих численных методов.
Продемонстрированная здесь возможность расчета НДС пластин с отверстиями различной толщины с помощью Ьр-МКНК имеет самостоятельное прикладное значение. Есть основания полагать, что развиваемый авторами МКНК может быть успешно применен для численного моделирования НДС многослойных пластин нерегулярных форм. При этом отличия значений характеристик НДС в различных теориях могут усиливаться при увеличении к, в том числе и восстановленных касательных напряжений, необходимых при изучении процессов разрушения композиционных материалов (см. [10, раздел 3.3.3]).
Возможность достижения высокой точности на грубых сетках и эффективность сочетания нескольких способов ускорения итерационного процесса решения получающихся в hp-МКНК СЛАУ открывает перспективы для проведения с его помощью оптимизационных расчетов НДС композитных элементов конструкций.
Работа выполнена в рамках государственного задания ИТПМ СО РАН. Авторы признательны А.Г. Горынину за полезное обсуждение результатов работы.
Литература
1. Ike, C.C. Mathematical Solutions for the Flexural Analysis of Mindlin's First Order Shear Deformable Circular Plates / C.C. Ike // Mathematical Models in Engineering. - 2018. -V. 4, № 2. - P. 50-72.
2. Тимошенко, С.П. Пластины и оболочки / С.П. Тимошенко, С. Войновский-Кригер. -М.: Физматгиз, 1966.
3. Reddy, J.N. Mechanics of Laminated Composite Plates and Shells: Theory and Analysis / J.N. Reddy. - Boca Raton; London; New York; Washington: CRC Press, 2004.
4. Голушко, С.К. Разработка и применение метода коллокаций и наименьших невязок к задачам механики анизотропных слоистых пластин / С.К. Голушко, С.В. Идимешев,
B.П. Шапеев // Вычислительные технологии. - 2014. - Т. 19, № 5. - С. 24-36.
5. Guo Chen. A Fast Finite Difference Method for Biharmonic Equations on Irregular Domains and Its Application to an Incompressible Stokes flow / Guo Chen, Zhilin Li, Ping Lin // Advances in Computational Mathematics. - 2008. - V. 29, № 2. - P. 113-133.
6. Ben-Artzi, M. An Embedded Compact Scheme for Biharmonic Problems in Irregular Domains / M. Ben-Artzi, J.-P. Croisille, D. Fishelov // Advanced Computing in Industrial Mathematics: 11th Annual Meeting of the Bulgarian Section of SIAM. - Cham: Springer, 2018. - Vol. 728. - P. 11-23.
7. Hailong Guo. A C0 Linear Finite Element Method for Biharmonic Problems / Hailong Guo, Zhimin Zhang, Qingsong Zou // Journal of Scientific Computing. - 2018. - V. 74, № 3. -P. 1397-1422.
8. Wenting Shao. Chebyshev Tau Meshless Method Based on the Integration-Differentiation for Biharmonic-Type Equations on Irregular Domain / Wenting Shao, Xionghua Wu, Suqin Chen // Engineering Analysis with Boundary Elements. - 2012. - V. 36, № 12. - P. 1787-1798.
9. Беляев, B.A. H-, p- и hp-варианты метода коллокации и наименьших квадратов для решения краевых задач для бигармонического уравнения в нерегулярных областях и их приложения / В.А. Беляев, Л.С. Брындин, С.К. Голушко, Б.В. Семисалов, В.П. Шапеев // Журнал вычислительной математики и математической физики. - 2022. - Т. 62, № 4. -
C. 531-552.
10. Идимешев, С.В. Модифицированный метод коллокаций и наименьших невязок и его приложение в механике многослойных композитных балок и пластин: дис. ... канд. физ.-мат. наук / Идимешев Семен Васильевич. - Новосибирск, 2016. - 179 с.
11. Garcia, O. hp-Clouds in Mindlin's Thick Plate Model / O. Garcia, E.A. Fancello, C.S. de Barcellos, C.A. Duarte // International Journal for Numerical Methods in Engineering. - 2000. - V. 47, № 8. - P. 1367-1522.
12. Baier-Saip, J.A. Shear Locking In One-Dimensional Finite Element Methods / J.A. Baier-Saip, P.A. Baier, A.R. de Faria, J.C. Oliveira, H. Baier // European Journal of Mechanics -A/Solids. - 2020. - V. 79. - Article ID: 103871. - 16 p.
13. Tsung-Hui Huang. A Stabilized Quasi and Bending Consistent Meshfree Galerkin Formulation for Reissner-Mindlin Plates / Tsung-Hui Huang, Yen-Ling Wei // Computational Mechanics. - 2022. - V. 70. - P. 1211-1239.
14. Слепцов, А.Г. Адаптивный проекционно-сеточный метод для эллиптических задач / А.Г. Слепцов, Ю.И. Шокин // Журнал вычислительной математики и математической физики. - 1997. - Т. 37, № 5. - С. 572-586.
15. Голушко, С.К. Метод коллокаций и наименьших невязок в приложении к задачам механики изотропных пластин / С.К. Голушко, С.В. Идимешев, В.П. Шапеев // Вычислительные технологии. - 2013. - Т. 18, № 6. - С. 31-43.
16. Vorozhtsov, E.V. On the Efficiency of Combining Different Methods for Acceleration of Iterations at the Solution of PDEs by the Method of Collocations and Least Residuals / E.V. Vorozhtsov, V.P. Shapeev // Applied Mathematics and Computation. - 2019. -V. 363. - Article ID: 124644. - 19 p.
17. Шапеев, В.П. hp-вариант метода коллокации и наименьших квадратов с интегральными коллокациями решения бигармонического уравнения / В.П. Шапеев, Л.С. Брындин, В.А. Беляев // Вестник Самарского государственного технического университета. Серия физико-математические науки. - 2022. - Т. 26, № 3. - С. 556-572.
18. Исаев, В.И. Варианты метода коллокаций и наименьших квадратов повышенной точности для численного решения уравнений Навье - Стокса / В.И. Исаев, В.П. Шапеев // Журнал вычислительной математики и математической физики. - 2010. - Т. 50, № 10. - С. 1758-1770.
Лука Сергеевич Брындин, младший научный сотрудник, Институт теоретической и прикладной механики им. С.А. Христиановича СО РАН (г. Новосибирск, Российская Федерация); младший научный сотрудник, центр НТИ по новым функциональным материалам, аспирант, кафедра математического моделирования, Новосибирский государственный университет (г. Новосибирск, Российская Федерация), [email protected].
Василий Алексеевич Беляев, кандидат физико-математических наук, инженер-исследователь, Институт теоретической и прикладной механики им. С.А. Христиановича СО РАН (г. Новосибирск, Российская Федерация), [email protected].
Василий Павлович Шапеев, доктор физико-математических наук, главный научный сотрудник, Институт теоретической и прикладной механики им. С.А. Христиановича СО РАН (г. Новосибирск, Российская Федерация), [email protected].
Поступила в редакцию 13 марта 2023 г.
MSC 35J40, 35J25, 35Q74 DOI: 10.14529/mmp230303
DEVELOPMENT AND VERIFICATION OF A SIMPLIFIED HP-VERSION OF THE LEAST-SQUARES COLLOCATION METHOD FOR IRREGULAR DOMAINS
L.S. Bryndin1'2, V.A. Belyaev1, V.P. Shapeev1
1Khristianovich Institute of Theoretical and Applied Mechanics SB RAS, Novosibirsk, Russian Federation
2Novosibirsk State University, Novosibirsk, Russian Federation E-mail: [email protected], [email protected], [email protected]
A new high-precision hp-version of the least-squares collocation method (hp-LSCM) for the numerical solution of elliptic problems in irregular domains is proposed, implemented,
and verified. We use boundary irregular cells (i-cells) cut off from the cells of a rectangular grid by a boundary domain and their external parts for writing the collocation and matching equations in constructing an approximate solution. A separate solution is not constructed in small and (or) elongated non-independent i-cells. The solution is continued from neighboring independent cells, in which the outer (and inner in a multiply-connected domain) part of the domain boundary contained in these non-independent i-cells is used to write the boundary conditions. This approach significantly simplifies the computer implementation of the developed hp-LSCM in comparison with the previous well-recommended version without losing its efficiency. We show reducing the overdetermination ratio of a system of linear algebraic equations in comparison with its values in the traditional versions of LSCM when solving a biharmonic equation. The results are compared with those of other papers with a demonstration of the advantages of the new technique. We present the results of bending calculations of annular plates of various thicknesses in the framework of the Kirchhoff-Love and Reissner-Mindlin theories using hp-LSCM without shear locking.
Keywords: least-squares collocation method; Kirchhoff-Love theory; Reissner-Mindlin theory; biharmonic equation; irregular domain.
References
1. Ike C.C. Mathematical Solutions for the Flexural Analysis of Mindlin's First Order Shear Deformable Circular Plates.Mathematical Models in Engineering, 2018, vol. 4, no. 2, pp. 50-72. DOI: 10.21595/mme.2018.19825
2. Timoshenko S.P., Woinowsky-Krieger S. Theory of Plates and Shells. New York, McGraw-Hill, 1959.
3. Reddy J.N. Mechanics of Laminated Composite Plates and Shells: Theory and Analysis. Boca Raton, London, New York, Washington, CRC Press, 2004. DOI: 10.1201/b12409
4. Golushko S.K., Idimeshev S.V., Shapeev V.P. Development and Application of Collocations and Least Residuals Method to the Solution of Problems in Mechanics of Anisotropic Laminated Plates. Computational Technologies, 2014, vol. 19, no. 5, pp. 24-36.
5. Guo Chen, Zhilin Li, Ping Lin. A Fast Finite Difference Method for Biharmonic Equations on Irregular Domains and its Application to an Incompressible Stokes Flow. Advances in Computational Mathematics, 2008, vol. 29, no. 2, pp. 113-133. DOI: 10.1007/s10444-007-9043-6
6. Ben-Artzi M., Croisille J.-P., Fishelov D. An Embedded Compact Scheme for Biharmonic Problems in Irregular Domains. Advanced Computing in Industrial Mathematics: 11th Annual Meeting of the Bulgarian Section of SIAM. Cham: Springer, 2018, vol. 728, pp. 11-23. DOI: 10.1007/978-3-319-65530-7_2
7. Hailong Guo, Zhimin Zhang, Qingsong Zou. A C0 Linear Finite Element Method for Biharmonic Problems. Journal of Scientific Computing, 2018, vol. 74, no. 3, pp. 1397-1422. DOI: 10.1007/s10915-017-0501-0
8. Wenting Shao, Xionghua Wu, Suqin Chen. Chebyshev tau Meshless Method Based on the Integration-Differentiation for Biharmonic-Type Equations on Irregular Domain. Engineering Analysis with Boundary Elements, 2012, vol. 36, no. 12, pp. 1787-1798. DOI: 10.1016/j.enganabound.2012.06.005
9. Belyaev V.A., Bryndin L.S., Golushko S.K., Semisalov B.V., Shapeev V.P. H-, p-, and hp-Versions of the Least-Squares Collocation Method for Solving Boundary Value Problems for Biharmonic Equation in Irregular Domains and Their Applications. Computational Mathematics and Mathematical Physics, 2022, vol. 62, no. 4, pp. 517-537. DOI: 10.1134/S0965542522040029
10. Idimeshev S.V. Modifitsirovannyy metod kollokatsiy i naimen'shikh nevyazok i ego prilozhenie v mekhanike mnogosloynykh kompozitnykh balok i plastin [Modified Method of Collocations and Least Residuals and Its Application in the Mechanics of Multilayer Composite Beams and Plates. PhD Thesis]. Novosibirsk, 2016. 179 p. (in Russian)
^.C. BpuHflHH, B.A. Be^HeB, B.n. OaneeB
11. Garcia, O. Fancello E.A., de Barcellos C.S., Duarte C.A. Hp-Clouds in Mindlin's Thick Plate Model. International Journal for Numerical Methods in Engineering, 2000, vol. 47, no. 8, pp. 1367-1522. DOI: 10.1002/(SICI)1097-0207(20000320)47:8<1381::AID-NME833>3.0.CO;2-9
12. Baier-Saip J.A., Baier P.A., de Faria A.R., Oliveira J.C., Baier H. Shear Locking in One-Dimensional Finite Element Methods. European Journal of Mechanics - A/Solids, 2020, vol. 79, article ID: 103871, 16 p. DOI: 10.1016/j.euromechsol.2019.103871
13. Tsung-Hui Huang, Yen-Ling Wei. A Stabilized Quasi and Bending Consistent Meshfree Galerkin Formulation for Reissner-Mindlin Plates. Computational Mechanics, 2022, vol. 70, pp. 1211-1239. DOI: 10.1007/s00466-022-02222-6
14. Sleptsov A.G., Shokin Yu.I. An Adaptive Grid-Projection Method for Elliptic Problems.
Computational Mathematics and Mathematical Physics, 1997, vol. 37, no. 5, pp. 558-571.
15. Golushko S.K., Idimeshev S.V., Shapeev V.P. Application of Collocations and Least Residuals Method to Problems of the Isotropic Plates Theory]. Vychislitel'nye tekhnologii [Computational Technologies], 2013, vol. 18, no. 6, pp. 31-43. (in Russian)
16. Vorozhtsov E.V., Shapeev V.P. On the Efficiency of Combining Different Methods for Acceleration of Iterations at the Solution of PDEs by the Method of Collocations and Least Residuals. Applied Mathematics and Computation, 2019, vol. 363, article ID: 124644, 19 p. DOI: 10.1016/j.amc.2019.124644
17. Shapeev V.P., Bryndin L.S., Belyaev V.A. The Hp-Version of the Least-Squares Collocation Method with Integral Collocation for Solving a Biharmonic Equation. Journal of Samara State Technical University. Series Physical and Mathematical Sciences, 2022, vol. 26, no. 3, pp. 556-572. DOI: 10.14498/vsgtu1936
18. Shapeev V.P., Isaev V.I. High-Accuracy Versions of the Collocations and Least Squares Method for the Numerical Solution of the Navier-Stokes Equations. Computational Mathematics and Mathematical Physics, 2010, vol. 50, no. 10, pp. 1670-1681. DOI: 10.1134/S0965542510100040
Received March 13, 2023
