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

Вестник КРАУНЦ. Физ.-мат. науки. 2011. № 2 (3). C. 18-32

Математическое моделирование Mathimatical simulation

УДК 534.121.2; 539.633 МОДИФИЦИРОВАННЫЙ МЕТОД ГРАНИЧНЫХ ЭЛЕМЕНТОВ ДЛЯ РЕШЕНИЯ ЗАДАЧ КОЛЕБАНИЯ ПЛАСТИН

Федотов В.П.

Институт машиноведения УрО РАН, 620049 г.Екатеринбург, ул. Комсомольская, 34

E-mail: fedotov@imach.uran.ru

Данная работа посвящена разработке алгоритмов решения волновых задач математической физики на основе метода граничных элементов (МГЭ). Основными преимуществами метода граничных элементов являются снижение на единицу размерности задачи и перенос дискретизации на границу исследуемой области, а также получение непрерывного решения внутри области. Вследствие этого сокращается объем расчетов и повышается точность решения.

Ключевые слова: метод граничных элементов, фундаментальное решение, граничные интегральные уравнения, напряжения, деформации.

(с) Федотов В.П., 2011

MSC 65N38

A MODIFIED BOUNDARY ELEMENT METHOD FOR SOLVING VIBRATION OF PLATES

Fedotov V.P.

Institute of engineering science Ural Branch RAS, 620049, Ekaterinburg, Komsomolskaya

st., 34

E-mail: fedotov@imach.uran.ru

This work is devoted to developing algorithms for the solution of wave problems in mathematical physics on the basis of the boundary element method (BEM). The main advantages of the boundary element method is reduced to one dimension of the problem and transfer sample to the boundary of the study area, as well as obtaining a continuous solution in the region. As a consequence, reduces the amount of calculation and improves the accuracy of the solution.

Key words: boundary element method, the fundamental solution, the oscillations, stress, strain.

© Fedotov V.P., 2011

Классификация численных методов

Исходная формулировка. Волновые задачи (к которым применим МГЭ [1]) можно свести к решению некоторого операторного уравнения

Ци) = Ь. (1)

Здесь u и Ь соответственно неизвестный и заданный вектора. Оператор Ь(-) в дальнейшем будет рассматриваться как дифференциальный оператор второго порядка, для уравнения колебаний он имеет вид:

д 2

Ц-) = 5(■) - д72 (■). (2)

Уравнение (1) можно записать в обобщенном смысле

J Ь(и(х))а (x)dx = J Ь(х)а (x)d х, (3)

о о

т.е. в (2) считается, что операторное уравнение (1) удовлетворяется не в каждой

точке рассматриваемой области О, а в целом по этой области с весом а(х). При

этом условия на границе должны выполнятся точно, что не всегда реально, либо приближенно:

ЗД = 8 х е Гь /(5(ц) -,у)О*(ауГь (4)

О(и) = g х е Г2,1 (О(и) - g)S*(ю)dГ2. (5)

Г2

где Г и Г2 - участки границы тела, в сумме дающие полную поверхность границы Г . Рассмотрим классификацию приближенных методов. Левую часть (2) называют внутренним произведением. В общем случае, дважды интегрируя левую часть равенства по частям, внутреннее произведение можно представить в виде

J Ь(и(х)) а(х)dх = J и(х)Ь*(а(х))dх + ^$*(а>)О(и) — S(а)О*(и))dx (6)

О О Г

где Г - внешняя граница области , О, S и О - дифференциальные операторы, обусловленные интегрированием по частям. Через S*(а) обозначены слагаемые, содержащие функцию а и получающиеся на начальном этапе интегрирования, S(u) содержит соответствующие слагаемые с функцией и . Оператор £*(■) называют сопряженным оператору Ь(^); оператор S(u) называется оператором существенных граничных условий, а оператор О(и) - оператором несущественных или естественных граничных условий. На границе области может быть задан любой тип граничных условий, но для единственности решения следует в некоторых точках задать существенные граничные условия, в частности для решения задачи колебаний необходимо хотя бы на части границы задать перемещения. Следующим шагом приближенного решения является переход от бесконечномерного уравнения (5) к конечномерному. Для этого искомую функцию ищут в виде

N

и(х) = £ а, ф,'(х), (7)

,=0

где <Рг(х) - известные функции, на которые накладываются условия, связанные с существованием решения. Система базисных функций <рг-(х), г = 0...^ должна быть линейно независимой, полной и должна удовлетворять граничным условиям. Выбор каждого метода зависит от выбора вида функций а(х). К исходной формулировке относятся метод коллокаций, метод наименьших квадратов, метод конечных разностей, метод Бубнова-Галеркина и др.

Слабая формулировка. Рассмотренные методы применимы для исходной формулировки (5), в которой участвует дифференциальный оператор второго порядка. Однако известно, что чем выше порядок производных в уравнениях, тем ниже точность приближенного решения. Это связано с погрешностями, возникающими при вычислении производной от приближенного решения. Для повышения точности обычно понижают порядок производной в исходном уравнении, т.е. переходят к слабой формулировке. Рассмотрим слабую формулировку для уравнения гиперболического типа. Используя формулы почленного дифференцирования и Гаусса-Остроградского, получим:

/ /(vu(x, т)у®(хт) - ~2д и(Хт) дадхт)) й Й(х) йт=

То а

Тр

= 1! ^и(х, т )а (х, т)йГ(х) й т — ^ ^ Т) ® (х, т) |^ ^(х) (8)

го Г Й

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

Обратная формулировка. Уравнение (7) дает слабую формулировку, которая снижает порядок производной на единицу в сравнении с исходной задачей, причем как искомая, так и весовая функции имеют первый порядок производной. Применение формул почленного дифференцирования и Гаусса-Остроградского к левой части (7) позволяет перейти к обратной формулировке, которая снижает порядок производных искомой функции еще на единицу за счет повышения порядка производной весовой функции. Выражение (3) преобразуется, для этого к первому слагаемому дважды применятся формула Гаусса-Остроградского - по границе Г(х)

Тр Тр

Jу"у2и(х, т)а(х, т)йЙ(х)йт = J !и(х, т)У2®(х, т)йЙ(х)йт+ (9)

То Г То Г

Тр Тр

+ J у"Уи(х, т )а (х, т)йГ(х) йт — J J и(х, т)У® (х, т)йГ(х) йт,

Т0 Г Т0 Г

второе слагаемое дважды интегрированием по частям:

Тр „ Тр

р [ 1 д2и(х, т) р г 1 д2а(х, т)

2 а (x, т) dQ(x) dт = u(x, т)^-—;—-— dQ(x) d т+

2 дт2 і і с2 дт2

c

to й to й

+ / “2дидГ1)®^ т) 1Тр йЙ(х) — /“2u(x, т) |Тр йЙ(х). (10)

Й Й Таким образом имеет место равенство

tF ' - 12

to й

V2u(x, т) —1 d u(x), т) ) а(x, т)dй^^т =

с2 дт2

V^(x, т) —— ^2’ ) u(x,tF)dQ(x)dт+ (11)

І д2а(x, т)

c

дт 2

Т0 Й

Тр Тр

+ J у"Уи(х, т)а (х, т) йГ(х) й т — J J и(х, т)Уа (х, т) йГ(х) й т+

Т0 Й Т0 Й

+ /“2 (т)и(х,т) — дидтт)а(х,т)) |т5 йЙ(х).

Й

Большим достоинством обратной формулировки является то, что производные здесь имеют место только для весовых функций, выбор которых зависит от исследователя. Методы граничных элементов, основанный на (11) получили развитие в работах [2]-[4], здесь в качестве весовой функции выбирается фундаментальное решение исходного уравнения. Для волнлвых уравнений фундаментальное решение - решение уравнения:

У2и*(х, () — 1д 2-,(| ^ Тр, т) = 3 (11£ — х||)8 (<р — т) (12)

представляет собой отклик среды на приложенное возмущение вида 8 - функции Дирака в точке £ в момент времени т. Фундаментальное решение (вследствие причинности) обладает свойством:

и*(£,х,Тр,т) = 0, / с(Тр — т) < |х — £| (13)

Применяя последнее свойство (13) фундаментального решения к выражению (11), а также учитывая, что верхний предел интегрирования обращается в ноль при Тр = т, поскольку и*(£,х,Тр,Тр) = ди*(£,х,Тр,Тр)/дТ = 0, для произвольной точки области можно записать соотношение:

Тр Тр

и(£, Тр ) = j j и(х, т ^*(£, х, Тр, т) йГ(х) й т — j j д(х, т )и*(£, х, Тр, т) йГ(х) й т+

Т0 Г Т0 Г

+ / С2 Г^u(^^tF,т)u*(xт) - ^т)u(^,x,tF,т)) |т=to d^x) (14)

й

где q(x,т) = Уи(х,т). В последнем соотношении интегрирование ведется по всей части границы, в то время как функции и(х,г) и q(x,г) заданы каждая только на своей части границы. Поскольку последнее уравнение справедливо на всей области Й ив том числе на ее границе Г функции и(х, г) и q(x, г) можно доопределить, взяв точку £ на границе и решив получившееся граничное интегральное уравнение:

где а(£) - функция угла, а(£) = 1/2 - для регулярной границы. Соотношение (14) и уравнение (1) справедливо для задач, сформулированных для областей произвольной размерности. Различия проявляются в виде фундаментального решения, определяемого уравнением (12).

ММГЭ для двумерного случая

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

здесь Г1 и Г2 = Г, Г1 П Г2 = 0 в совокупности дают границу Г области Й.

Фундаментальное решение. Фундаментальное решение в двумерном случае -отклик среды на приложенное в точке £ и момент времени т возмущение вида 8 -функции представляет собой круговую волну с центром в точке £

tF

tF

to Г

to Г

й

дudx т) u(£, x, tF, т)^ |т=to da(x), (15)

(16)

w(x, to) = Ф (x), Wt (x, to) = Y (x),

(17)

с граничными условиями:

w(x, t) = w(x, t) x є Г1,

(18)

d (x, t)

—г— = q(x, t) x є Г2,

д n

(19)

M*(£, x, tF, т)

cH (c (tF — т) — r)

c(tF — т) — r = o,

(20)

2 n \J c2 (tF — т )2 — r2 ’

где Н(х) — ступенчатая функция Хевисайда, г = |х — £ | — расстояние от точки влияния х до точки £. Производная по нормали к границе Г от функции и-(£,х,Г',т)

д-(£, х, Г', т) = ^П- = —С1'Н (С (Г'—-—Г)^у2 Пг, с(Г' — т) — г = 0, (21)

дП - С2 ^2_ г^3/2

2п (^с2 (Г' — т)2 — г2

где пг = дг/дп производная длины радиус-вектора по нормали к границе, а производная по времени от фундаментального решения

д и- (£, х, Г', т) с2 (Г' — т) Н (с (Г' — т) — г)

дт ( 2 ч-з/Г, с(г'— т) — г =0 (22)

2п (с2 (Г' — т)2 — г2^

Граничное интегральное раннение. Граничное интегральное соотношение (15), с учетом вида фундаментального решения (20) и его производной по нормали на границе (21), принимает вид:

и(£, Г')= I' (и(х, т)—гН(с(г'—^—г)__Пг^Г(х)^т—

' Г 2п(с2(^ — т)2 — г2)3/2

— // Ув(х, т) Н (с(>' ~ т) ~ г) ^ Г(х)^ т + (23)

I Г ( , ) 2п^с2(<Г — т)2 — г2 ( ) ' '

+ 1 гафт)|т=,о Н(с(»—<0) — г) ^п(х) + А /с2(<г—<0)Н(с(»—,0)-г)Лй(х).

С й дт 2п V'""С:(/. — Го)2 — г2 С й 2п (с2((^ — г0)2 — г2) /

Первый и последний интегралы, входящие в это соотношение, по отдельности не сходятся, поскольку имеют неинтегрируемые в смысле главного значения особенности в окрестности точек с(Г' — т) — г = 0, , однако их сумма конечна, поскольку сходятся остальные интегралы, входящие в соотношение. В отдельных случаях граничное интегральное соотношение преобразуется таким образом, что неинтегриру-емые особенности сокращаются. Начальные условия заданы в виде гармонических функций. В этом случае с помощью второго тождества Грина [1] интеграл по области можно преобразовать в эквивалентные граничные интегралы. Для случая, когда П - гармоническая функция в области а , т.е. имеет место равенство Ап = 0 , можно записать вторую формулу Грина в виде:

I (пАЦ, + Ц)Ап)^ = | ГпддиП0 — цдПЛ (24)

а г

Кроме того, если функция Ц0 подобрана так, что АЦ0 = и0, то

/П“0<® = /(-£ — Ц0 Щ ) ^Г (25)

а г

Одну из таких функций можно найти из соотношения

^0(г) = У1 ^ гю(г)^^ ^г (26)

Простой подстановкой можно убедиться, что:

дм*(£, x, tF, т) { c (tF —10) H(c (tF —10) — r)

эГ k=t0 = — „ / _ 2\3/2

2п (c2 (tf — to)2 — r2j

72^ cH(c (tF — t0) — r) l^2 (c (tF — t0) — RF A + _ ^2Г

2 п I r /2 п

где = а/с2 (Г' — Г0)2 — г2. Производная по нормали к границе от функции и- равна: ди- | с2 (Г' — ?0) Н(с (Г' — ?0) — г) 1

r

дп 1 - I 2 и 7^2 2 2пг

' 2 п гл/с2 (Г' — Г0) — г2

Применяя второе тождество Грина к последнему из интегралов правой части выражения (22), получаем:

1 (г (х Г ) с2(г' — Г0)Н(с(г' — Г0) — г) л0(х)

С2 и(х , Г0)-----(-------------------_2- ^а(х) =

а 2п (с2 (Г' — Г0) — г2)

[ ( .4 с2(г' — Г0)Н(с(г' — Г0) — г) , 1 ) , (27)

= —У и(х,,0) I . /2, ^~=T~ + 2пг ) МГ(х)+ (27)

2п r yc2 (tf - /ні2 - r2

+.1 I дЦ0 СН(с (Г' — Гр) — г) / 2 (с (Г' — Г0) — А + 1п(гЛ ^Г

с2У дп \ 2п \ г / 2п /

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

д | (Г' — т) Н (с (Г' — т) — г) ) гН (с (Г' — т) — г)

— /-2 (^_ т)2- г2 "г = 2 п (с2 ^2 г2Ч3/2 Пг,

2пгус2 (Г' — т) — г2 ) 2п ^с2 (Г' — т)2 —

тогда в первом слагаемом выражения (23) интегрирование по времени проводим по частям:

гН(с(г' — т) — г) , Л [ , Л (г' — Г0) Н(с(г' — Г0) — г)

to Г

rH(c(tF — т) — r) Л, f , Л (tF —ч-r *oy -у Л

m(x, т)----------------------3/2 nrd r(x)d т = m(x, t0)------- nrd r(x) —

2п (c2 (tF — т)2 — r2j г 2пrye2 (tf — to)2 — r2

f f Sт) (/f - т)И(c(/f - т) - r)nrdr(x)йт

,0 Г дт 2пг\!с2 (<' — т)2 — г2

Здесь учтено, что верхний предел интегрирования обращается в ноль при Г' = т, и проаедена смена местами пределов интегрирования во втором слагаемом правой

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

ч [ С д и(х, т)(?' — т) Н (с(Г' — т) — г) тч-,/ \

и(£, Г') = — —^Пг^г(х)^ т—

( Г дт 2пг^/с2 (* — т)2 — г2

— / / У“(х, т) 2Н— т) — г)=т^Г(х) Лт+

У Г 2п^/с2 (Г' — т)2 — г2

+ 1 г д^ Н(С (,' — >0) — г) ^Й(Х)+

С2 а дт |т 2пV2 (Г' — <0)2 — г2 ( )

+ І г Гдm(x, to) cH(c(tf — to) — r) l/2(c(tf — to) — Rf ))\ + ln(r) + c2J d n 2п П r )+ 2 п

Г

/n

m(x, to )^-^d r(x).

2 п r

dr(x) +

Пг

2пг

Все интегралы, стоящие в правой части последнего выражения, интегрируемы в обобщенном смысле. Кроме того, благодаря свойству гармонических функций, интегрирование функции, определяющей начальное отклонение мембраны, по области а удалось заменить эквивалентными интегралами по границе Г. Для того, чтобы избежать интегрирования по области, интеграл от функции, характеризующей начальные скорости точек мембраны, можно также преобразовать в эквивалентные граничные интегралы для случая, когда ди(х, т)/дт|т=Г0 - гармоническая функция. Простой подстановкой можно убедиться, что:

и-(£, х, Г', Г0) = У2Ц-

и - = Н (С (гг — ^ — г) (с — ) ( 2 (С (Г' — Г0) - *')) — 0) 1п(г)

2 п V V гс2 (Г' — Г0)2 } ' ) 2 п 47

Производная по нормали к границе:

ди*_ Н(с (Г' — Г0) — г) ( г | С (Г' — *0)^ П | С (Г' — Г0) П

+ ) Пг + _ Пг

д п 2 п — с (Г' — Г0) г ) 2 п г

Применяя второе тождество Грина к интегралу, характеризующему влияние начальных скоростей мембраны, получаем:

1 г йй(Х) =

с2 а дт |т Г02п0)2—3 ()

1 [ ди(х, тV ди- 1 [ д (ди(х, т)

С2 У дт 0 дn С2 J дn V дт

ГГ

d^ |т=to U*(^,x,/f,to)dr(x)

Учитывая начальные условия (16), граничное интегральное соотношение, в случае, когда начальные условия заданы в виде гармонических функций имеет вид:

Таким образом, перемещения выражены только через функции, заданные на границе, следовательно, для их нахождения необходимо разбиение на элементы не всей области П, а только ее границы Г. Размерность задачи снижена на единицу.

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

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

Выберем теперь элемент границы - отрезок. Относительно него введем специальную систему координат таким образом, что ее начало O совпадает с началом отрезка, ось ОХ направлена вдоль направления, задаваемого отрезком, ось ОГ образована поворотом оси OX на угол 2 против часовой стрелки. Такой элемент в дальнейшем договоримся называть базовым, а систему координат - системой координат,

tF

дu (x, т) (tF — т) H (c (tF — т) — r)

дт

nrd Г^) d^

to Г

tF

J У Vu(x т)

(28)

to Г

+cL у ^(x)^ d Г^) + cL у U *(£, x, tF, to)d^x)+

ГГ

Г

Г

связанной с базовым элементом. Направлением нормали будем считать направление, противоположное оси ОГ. При выводе граничного интегрального соотношения подразумевалось, что нормаль направлена «наружу» области. Таким образом, при рассмотрении области, ограниченной кусочно-линейной границей, считаем, что задача решается внутри области, если обход границы (ориентация отрезков) ведется против часовой стрелки; и вне области, если обход границы ведется по часовой стрелке.

В данной работе приведена реализация метода, с линейной по координате интерполяцией внутри базового элемента значения деформаций и перемещений. Сделаны следующие предположения. Пусть значения функции и(х, т) аппроксимированы внутри пространства [ту—1, ту] х Гг- полиномом первой степени по координате и времени, построенным по четырем точкам: (хг-_1,т/-_1), (хг-_1,ту), (хг,ту_1), (хг,ту).

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

В случае, когда два соседних элемента границы расположены под углом, отличным от п, они имеют различные нормали и, соответственно, в окрестности смежной

д и

с ними точки различные производные по нормали —.

д п

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

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

Для согласования значений напряжений и деформаций на соседних граничных элементах необходимо вычислять номера соседних элементов в специальной нумерации. Для вычисления номеров этих элементов введены функции /1 : N ^ N - возвращает номер в прямой нумерации по номеру в специальной и /2 : N ^ N - возвращает номер в специальной нумерации по номеру в прямой. Таким образом, в специальной нумерации соседний правый элемент г-го элемента границы будет иметь номер /2(Ж0 + 1), а соседний левый /2(/1(г) — 1). Для краткости соседний правый элемент в специальной нумерации будем обозначать индексом г+, а соседний левый индексом г—.

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

такой узловой точке необходимо знать одно значение для и и два значения для д, которые будем обозначать д^1 ) - значение для текущего граничного элемента в на-

(3)

чальной точке и д;:_і - значение для предыдущего граничного элемента в конечной

точке. Аналогично, для функции и будем обозначать и(.1) и и(3). , которые очевидно совпадают. Опуская промежуточные выкладки, получим для системы (28) следующую конечно-разностную схему

1 N , ч , ч 1 ^+М , ч , ч

V У"д(р,0)^(р,0) + V V и(р,1)^(р,0) —

II дгк г1,шгкк + I I игк г2,шгкк —

р—0 г—1 р—0 г—N+1

— В(^И) Г; 'к, '0) с(£ш)и(£ш5 'к) 5 т — 15.. ^5

1 N / \ / ч 1 N+М

с(«.)и(6.5*) + II+1 I и^™ —

р—0 г—1 р—0 N +1

— В(£ш, Г, , ?0), . — N + 1,..^ + М,

. — ^1(1,0)(«5 Г'5 5 О-ь О')5 ^5 Г'5 'ь О-ь О)5

1 N 1 N+ra

3(5.5 Г 'к 5 '0) — - I I иГ^ - I I -

р—0 г+1 Р—0 г—N+1

к-1 N+M / 1 , ч 1 .А

11 I ^25.^ +1 *0І +

г—1 і—1 \р—0 р—0 /

+ _! У у (Vr)F(r) + W(r)F(r) + ©WF(r) + ©(r)F(r)

+ c2 І І [ті F3,mik + Ги,г M,raik + ^Пг F5,mik + F6,mik

1 1

С2

с г—1 г—1

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

Примеры

Пример 1. Мембрана имеет форму прямоугольника АВСО. Размеры мембраны АВ = СО = Ь, ВС = ОА = 1. Стороны мембраны закреплены, причем в интервале времени ? = [0, п] движутся по закону и = — 8т(?)/10, и далее неподвижно фиксируются в нулевом положении. Начальные отклонения и скорости точек мембраны нулевые.

Размеры мембраны, выбранные для решения задачи Ь = 1 = 5 ед, скорость звука принята равной 1. Задача решалась до 20 с. Рассмотрены три варианта решения с

различными шагами по времени дЛ и границе ^й. Первый вариант дЛ = 1 с, = 2 ед, второй вариант Л = 0.5 с, = 1 ед, третий вариант Л = 0.125 с, = 0.25 ед. Колебания центральной точки х = 2.5;у = 2.5 изображены на Рис. 1.

Рис. 1. Отклонение центральной точки квадратной мембраны. Размеры мембраны 5 на 5 ед. Стороны мембраны в интервале времени г = [0, п] движутся по и = — зт(г)/10, и далее фиксируются в нулевом положении.

Из рисунка видно, что хорошая сходимость обеспечивается уже во втором варианте решения ( Л = 0,5 с, = 1 ). Расчеты, проведенные классическим МГЭ, совпали с результатами ММГЭ. Время решения задачи с разбиением Л = 0,5 с,

= 1 для МГЭ составило 46,68 с, для ММГЭ 0,78 с. Таким образом, при решении задачи ММГЭ за счет однократного вычисления интеграла по базовому элементу достигнуто существенное ускорение счета.

Сравнение времени и точности решения с классическим МГЭ приведено в таблице.

Метод решения Невязка Время решения, c

МГЭ (Л = 0,5, ^ = 1,0) 0,05 4,687

ММГЭ (Л = 0,5,^ = 1,0) 0,04 0,781

МГЭ (Л = 0,125,^ = 0,25) 0,0028 346,158

ММГЭ (Л = 0,125,^ = 0,25) 0,0021 3,126

Несколько кадров решения показаны на Рис.2.

Пример 2. Мембрана имеет форму прямоугольника АВСО, стороны мембраны ВС и ОА не закреплены и свободно движутся, а стороны мембраны АВ и СО закреплены, причем сторона СО постоянно закреплена в нулевом положении, а сторона АВ совершает одно полуколебание по закону и = 8т(г)/10 в интервале времени г е [0, п] и далее фиксируется в нулевом положении.

При этом возникает волна, распространяющаяся от АВ. Задача обладает симметрией и сводится к одномерному случаю, рассмотренному в разделе 2.2, таким образом, имеется точное решение, с которым и будем сравнивать решение, полученное ММГЭ.

Размеры мембраны, выбранные для решения задачи АВ = СО = 10 ед, ВС = ОА = 5 ед, скорость звука принята равной с = 1 . Задача решалась до 20 с. Рассмотрено

Рис. 2. Отклонения мембраны.

два варианта решения с различными шагами по времени Л и границе ^й. Первый вариант Л = 0.5 с, = 1 ед, второй вариант Л = 0.25 с, = 0.5 ед. Колебания точек (х = 1.25;у = 0.5) и (х = 1.25;у = 0.5) для первого варианта решения изображены на Рис.3.

Рис. 3. Отклонения точек мембраны.

Колебания (х = 1.25;у = 0.5) и (х = 1.25;у = 0.5) для второго варианта решения изображены на Рис.4. Непрерывной линией показано точное решение.

Рис. 4. Отклонения точек мембраны.

Рис. 5. Отклонения мембраны. Справа точное решение.

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

Заключение

1. В отличие от классического МГЭ в ММГЭ интегрирование проведено аналитически, причем не по каждому элементу границы, а по единожды выбранному базовому элементу, помещенному в начало координат. Благодаря этому интегрирование по элементам границы сводится к подстановке координат точки в компактные формулы, выраженные через элементарные функции.

2. В случае, когда шаги по времени являются постоянными, матрица разрешающей системы не меняется, и решение СЛАУ сводится к простой операции умножения матрицы на вектор.

3. Один раз решив задачу определения граничных значений, мы можем сколько угодно раз менять сетку внутренних точек и пересчитывать необходимые значения величин в них без дополнительного составления и решения СЛАУ.

4. На всех этапах решения задачи произведено полное распараллеливание счета, лежащее в основе ММГЭ. На первом этапе формируется система линейных уравнений, при этом распараллеливание счета является абсолютным, так как вычисление коэффициентов разрешающей системы осуществляется подстановкой координат точек влияния в элементарные функции независимо друг от друга. На втором этапе решается СЛАУ, при этом матрица разрешающей системы остается неизменной на всех шагах по времени. Поэтому собственно решение СЛАУ происходит только на первом шаге, а на последующих шагах решение системы сводится к быстрой операции умножения матрицы на вектор. Такая операция эффективно распараллеливается. Первые два этапа повторяются на каждом шаге по времени. Наконец, на третьем этапе вычисляются перемещения точек мембраны внутри области путем подстановки их координат в элементарные функции. Расчет перемещений точек, как и на первом этапе, производится независимо, что позволяет абсолютно распараллелить вычисления.

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

Литература

1. Тихонов А. Н., Самарский А. А. Уравнения математической физики. М.: Наука, 1977. 728 c.

2. Бреббия К., Теллес Ж., Вроубел Л. Методы граничных элементов M.: Мир, 1987. 524 c.

3. Mansur W.J., Brebbia C.A. Formulation of the boundary element method for transient problems governed by the scalar wave equation // Appl. Math. Modelling. 1982. Vol. 6, iss. 4. P. 307-311.

4. Федотов В. П., Спевак Л. Ф. Решение связных диффузионно-деформационных задач на основе алгоритмов параллельного действия. Екатеринбург: УрО РАН, 2007. 191 c.

Поступила в редакцию / Original article submitted: 03.10.2011