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

МЕХАНИКА

УДК 539.3

РЕШЕНИЕ ТРЕХМЕРНЫХ ЗАДАЧ ДИНАМИЧЕСКОЙ ТЕОРИИ ПОРОУПРУГОСТИ МЕТОДОМ ГРАНИЧНЫХ ЭЛЕМЕНТОВ С ПРИМЕНЕНИЕМ ПАРАЛЛЕЛЬНЫХ ВЫЧИСЛЕНИЙ

© 2011 г. Л.А. Игумнов, И.С. Карелин

НИИ механики Нижегородского госуниверситета им. Н.И. Лобачевского

[email protected]

Поступила вледакцою 03.03.2011

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

Ключевые слова: среда Био, метод граничных элементов, трехмерные задачи динамики, параллельные вычисления.

Введение

Уникальные свойства пористых материалов послужили их широкому распространению в машиностроении, строительстве и т.п. Для исследования распространения волн в пористых средах наиболее часто используется модель пороупругой среды, предложенная Био [1-4]. К подобным исследованиям может применяться метод граничных элементов (МГЭ) [5, 6]. Привлечение гранично-элементных технологий к численному моделированию динамики пороуп-ругих тел находится в фазе активного становления. Построение современных численных схем МГЭ невозможно без использования параллельных вычислений.

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

Математическая модель и граничноинтегральные уравнения

Система дифференциальных уравнений в преобразованиях Лапласа (параметр 5) для смещения и- и порового давления р имеет следующий вид [6]:

Ои<Л + ^К + 3 О)иММ -

-(а - Р) р- 5 2(р - Рр /)и,- = - ^,

— рм - ^ Р - (а - Р)5иу = -а ,

5р / К

Р = ф2S + S2k(ра + фр/ ^ где О, К - константы материала, К - константа модели материала Био, ф - пористость, k - проницаемость, а - эффективный коэффициент напряжений, р, ра, р/ - пористости скелета, присоединенной массы и жидкой среды, Р{,а -плотности источников.

Интегральное представление прямого подхода имеет следующий вид [6, 7]:

йГ -

1 <;з I = [ - 1 - 1 1

1 1 J Г и/ 1 - 1 1

- &" 1 1

Т/ / _ _р _

йГ.

Ядра допускают выделение особенностей Р’=0(г0),

и/ =0(г0),

-{г/м + 5„ (3 - 4У)}1 + 0(г0),

4пР г

<2

1 + V

-{а(1 - - п-) -

' 8п(1 - V)

- 2Р(1 - 'ОО^- +п' )}1 + 0(г 0),

/ Р/Р Ь пч1 - 2v

Т1 =~^ \ (а - Р)----г Г' +

1 - V

+ п

8лР

а + Р(1- 2у) 1 ^ + 0(г 0),

т

1 - V

, _ - [(1 - 2^5„ + 3Г/,' Ь +

8п(1 - v)r2

(1 - 2v)(r'.п,. - ли.) „

+ V------, 2^ о(г0),

8п(1 - v)r

л Г

О1 = --д+0(г0).

Итоговая система гранично-интегральных уравнений (ГИУ) имеет вид:

dГ =

Ч-(у) 0 1 и , (я, х)

_ 0 с(у)_ _Р(S, х) _

~Т8 (8, у, X) - °8(Х y, х) и., (я, х)

т/(8 У, х) 1 ,х) ,у, - _ Р(S, х) _

Дл>1

dГ.

/ (0)

2п

Ек+\ Ек (с°5 (га*+1?)-

Д 7,

где

- cos (шкг))+ °к+Д Ок (ап (ш*+1?)-Д к

- sin (га к?))]

дк = гак+1- гак> рк = Re[ 1 (а +,тк)]

1 / , т = Т 8 - °"

1 / ‘С. 1 / - о1 _

Базовый процесс ГЭ-дискретизации состоит в разбиении поверхности сЮ на NE граничных элементов Ее (1 < е < NE). Граничные элементы

- четырехугольные восьмиузловые биквадра-тичные (рис 1). На каждом элементе вводится локальная система координат £ = (£ ,£2)е

е [-1,1]2. Глобальные координаты произвольной точки элемента Nk выражаются через локальные с помощью функций формы. Для преобразования локальных координат точки в глобальные используется следующее выражение:

(££) = Е N (££)

в( к‘ \, = 1,2,3,

£ е Д е ,(2)

и 8 (я, у, х) - Р' (я, у, х) (s, х)

и( (я, у, х) - Р/(я, у, x)\\_q(s, х)_ решения проблемы обращения интегрального преобразования Лапласа применим вариант метода Дурбина [8]:

^+, - ^)Дк

гдеР(к, I) - глобальный номер узла, имеющего в к-м элементе локальный номер I, N (£) -функции формы. В качестве функций формы выбраны квадратичные полиномы интерполяции.

Неизвестные граничные поля (и, ?) интег-

к / к \

рируются через узловые значения и = и(2 ) и ~к = ~ (хк) в интерполяционных узлах 2к. Множество интерполяционных узлов отличается от множества геометрических узлов, а множество интерполяционных функций не совпадает с множеством функций формы. Рассмотрим случай, называемый согласованным интерполированием (Р.В. Гольдштейн (1978)), где для аппроксимации граничных перемещений применим билинейные элементы, а для аппроксимации поверхностных сил - постоянные элементы. При этом для расчетного значения параметра р будем иметь следующие выражения граничных перемещений и поверхностных сил внутри элемента Sk:

, (у) = 2 X (£)и

Х(к ,1)

Ок =1т[ 1 (а + /х к)].

Гранично-элементная дискретизация

Чтобы ввести гранично-элементную (ГЭ) дискретизацию, рассмотрим регуляризованное уравнение:

апик(х) + [ »~к(x, У, 8)к-(у) - (х)] -

сЮ (1)

- Ол (x, y, 5)Х(у)}йх = o, (х е ащ

где

й/1 ~ \т? - о8"

О = . , Т = Т8 ,Т = [?^]т.

+

+

и

Рис. 2

Рис. 3

~ (у) = ^ і = 1,2,3; у є Sk.

Здесь Я1 (£) - функции формы для линейного четырехугольного элемента.

Для получения дискретного аналога ГИУ применим метод коллокации. В качестве узлов

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

1 _ N 4

ит+ ЕЕ а

т,к,1

2

их(к,/) =

к=1 =1

N

(3)

Е в.

і ч

к=1

Л N 4

ЬОЛ и т +

>т + ЕЕ і

к=1 1=1

N

= Е втк тк

і ]

тЫ их( к,/) =

(4)

где N - число элементов границы.

Уравнения (3) записаны в узлах аппроксимации обобщенных перемещений, уравнения (4) записаны в узлах аппроксимации обобщенных усилий:

1 1

' = Ц[Д' (£)-(хт,ук Ш р) --1 -1

- V ),тТ0 (хт,ук т р)У 1 1

В? = 11 и8 (хт,ук Ш рУк О&УШг .

-1 -1

Необходимо отметить, что коэффициенты дискретных аналогов (ДА) имеют особенность типа 1/г и 1/г2. Это определяет специфику вычислительного процесса. При вычислении интегралов по поверхности рассматриваются два

варианта. Первый - точка хт не принадлежит элементу интегрирования. В этом случае интегрирование по элементу сведено к повторному интегрированию по локальным координатам £1 и £2. По каждой из координат используются квадратурные формулы Гаусса. Второй вариант

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

Для каждого уравнения системы (3)-(4), зато

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

Рассмотрим построение параллельного алгоритма вычисления коэффициентов для уравнений (3). Параллельный алгоритм строится на базе последовательного [9]. Блок-схема последовательного алгоритма представлена на рис. 2.

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

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

8

к=1

У О У О

/=Зм Рис. 4 409

Рис. 5

At=0.0001 с

t,c

Рис. 6

t,c Рис. 7

узлов аппроксимации обобщенных смещений. Для вычисления коэффициентов ДА уравнений (4) применяется аналогичный алгоритм. Его степень параллелизма равна числу узлов аппроксимации обобщенных усилий.

Параллельный алгоритм реализован на языке Фортран с использованием библиотеки параллельных вычислений ОрепМР [10] для реализации потоков исполнения. На этапе работы программы определяется количество ядер процессора и запускает такое же число потоков исполнения для достижения наибольшей производительности на данной ЭВМ.

Численные результаты

Рассмотрена задача, изображенная на рис. 4, со следующими параметрами материала:

E = 1.44 • 1010 H / м 2;v = 0;р = 2458 кг / м3; р f = 1000 кг / м 3;ф = 0.19;

R = 4.7 • 108 H / м2; а = 0.86; k = 1.9 • 10 10 м4/Нс.

Нагрузка, действующая на тело, t =

= -1 Н/м2.

Граничные условия в области Лапласа имеют вид: йу (у = 0) = 0, qу (у = 0) = 0, оy (y = 1) =

= -1, p(y = l) = 0.

Гранично-элементная сетка, изображенная на рис. 5, состоит из 504 элементов. Результаты расчетов перемещений и напряжений приведены на рис. 6, 7 соответственно.

Вычисления производились на двухъядерном процессоре Intel Core 2 Duo CPU E8400. В таб-

Таблица

Затраченное время, Tp

Алгоритм Последовательный 56 112

Параллельный 21 4З

Ускорение алгоритма, S 2,67 2,60

Количество частот 5 10

лице представлено время, потраченное программами на расчет решения для 5 и 10 частот. Для вычисления ускорения алгоритма Sp применялась формула:

где 7], Tp - время выполнения программы на

одном и на p процессорах соответственно. В сравнении с последовательным алгоритмом реализация параллельного дает значительный прирост производительности.

Список литературы

1. Френкель Я.И. К теории сейсмических и сейсмоэлектрических явлений во влажной почве // Изв. АН СССР. Сер. географ. и геоф. 1944. Т. 8. № 4. С. 65-78.

2. Biot M. Theory of propagation of elastic waves in a fluid-saturated porous solid. I. Low-frequency range // J. Acoust. Soc. Am. 1956. V. 28. № 2. P. 168-178.

3. Biot M. Theory of propagation of elastic waves in a fluid-saturated porous solid. II. Higher-frequency range II J. Acoust. Soc. Am. 1956. V. 28. № 2. P. 179— l9l.

4. Николаевский В.Н. Геомеханика и флюидо-динамика с приложениями к проблемам газовых и нефтяных пластов. М.: Недра, 1996. 447 с.

5. Дунин С.З., Михайлов Д.Н., Николаевский В.Н. Продольные волны в частично насыщенных пористых средах. Влияние газовых пузырьков II ПММ. 2006. Т. 70. Вып. 2. С. 282-294.

6. Schanz M. Wave propogation in viscoelastic and poroelastic continua. Berlin: Springer, 2001. 170 p.

7. Аменицкий А.В. Развитие метода граничных элементов для численного моделирования динамики трехмерных однородных пороупругих тел: Автореф. дис.... канд. ф.-м. наук: 0l.02.04. Н.Новгород, 2010. 20 с.

В. Баженов В.Г., Белов А.А., Игумнов Л.А. Гранично-элементное моделирование динамики кусочно-однородных сред и конструкций: Учебное пособие. Н.Новгород: Изд-во Нижегородского госунивер-ситета, 2009. 1В0 с.

9. Breshears C. The Art of Concurrency. O’Reilly Media, Inc., 2009. З02 p.

10. Chapman B., Jost G., Van der Pas R. Using OpenMP. Portable Shared Memory Parallel Programming. The MIT Press Cambridge, Massachusetts, London, England. 200В. З5З p.

SOLVING 3D PROBLEMS OF DYNAMIC POROELASTIC THEORY BY BOUNDARY ELEMENTS

METHOD USING PARALLEL COMPUTING

L.A. Igumnov, I. S. Karelin

Direct boundary element method for solving 3D dynamic poroelastic problems is considered. Parallel computation method is used to increase efficiency of numeric modelling. Results of computation experiments are presented.

Keywords: BEM, Biot porous media, boundary elements, boundary integral equations, parallel technology.