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

УДК 531/534

В.П. Федотов, Л.Ф. Спевак, В.Б. Трухин, В.В. Привалова, Т.Д. Думшева, Е.С. Зенкова

ИССЛЕДОВАНИЕ СХОДИМОСТИ ЧИСЛЕННО-АНАЛИТИЧЕСКОГО МЕТОДА РЕШЕНИЯ ЗАДАЧ УПРУГОСТИ, ТЕПЛОПРОВОДНОСТИ И ДИФФУЗИИ

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

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

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

1. Алгоритм решения упругой задачи. Описание численно-аналитического метода решения, основанного на методе граничных элементов, начнем с задачи теории упругости. Рассмотрим плоскую упругую область произвольной формы, на границе которой заданы некоторые нагрузки или перемещения. Задачей является нахождение вектора перемещений и1, тензора деформаций в - и тензора напряжений с-, которые удовлетворяют в области системе уравнений:

сш =Ь-; (1)

еИ =(*—■ + иг,} V2; (2)

С- = 2теу + ^(Вп + 822 К- (3)

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

на поверхности Бу : спп1 + ст12п2 = /1 = /*, ст21 п1 + с22п2 = /2 = /2*;

на поверхности Би: и1 = и1 , и2 = и2 . (4)

Здесь ЬI - известные объемные силы, с-г- = 5с1 -1дх1 + дс2дх2 , здесь и далее по повторяющемуся индексу производится суммирование от 1 до 2, и{ - = диг /дх- , т - модуль упругости при сдвиге, V - коэффициент Пуассона, 5- - единичный тензор, п - вектор нормали к поверхности. Обратная к (3) зависимость будет иметь вид

8« = 2г(С"'“ +С 22 Ч . (5)

Согласно методу граничных элементов [1, 2] решение системы (1) - (4) можно выразить через поверхностные перемещения ui (х) и поверхностные напряжения (х) . Для вектора перемещений в любой внутренней точке можно записать

иг (х) = | к (Х х)/] (х) - (Х х(х^ + | к (Х х)/- (х) - (Х х)и* (х)К (6)

и

Здесь х е Б - граничная точка области, X - внутренняя точка области, звездочкой (*) обозначены известные из граничных условий для группы поверхностей Би перемещения и * (х) и для группы поверхностей Б/ поверхностные напряжения /*(х). Функции и* (X, х) (тензор Грина) и /*(£>, х) для двумерной задачи имеют вид:

и* (Х х) = о_(, 1п)..[(3 - 4п)1п (г )5У -D irD / ],

8л(і - у)т - і

^(Х) = 4р(-у) {[- 2п) + 2ЛгГА}Г] - (і - 2п)'ГП] - Л]Шг ^:

где г = г(х,X) - расстояние между х и X , Ліг = .

дХ:

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

2 иг (х0 ) = | ^ (х0 , х)/] (х) - /Ц (х0 , хК (х)]& + | [и* (х0 , х)/1 (х) - /Ц (х0 , хК (х)]& (7)

2

и используется для определения неизвестных поверхностных перемещений и - (х) и поверхностных напряжений /- (х). Таким образом, можно интерпретировать изучаемый процесс как

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

Разбив границу области на элементы (рис. 1), мы приходим к численной схеме, в основе

которой лежит участок граничного элемента (отрезок прямой либо дуга окружности) и произвольная точка, на которую оказывают влияние напряжения и перемещения, действующие на этом элементе. Если сделать разбиение поверхности деформируемой области достаточно мелким, то все поверхностные перемещения и напряжения, включая неизвестные, можно считать постоянными в рамках отрезков и помещенными в середину элемента. В этом приближении, записав уравнение (7) для середины каждого граничного элемента, получим следующую систему уравнений:

(хі; х2)

Р и с. і. Область произвольной геометрии и влияние элемента на точку

■ +

+

2 и() = Е /Р |и* х)*-Е и ;а) I fУ(x(P), х)/5-

2 а=1 с (а) а=1 с (а)

С/ С/

N +М / ч N +М / ч

Е /у°° |и* (х(Р), х)* - £и(“> I/г*(х(р), х)& , г = 1,2, р = 1..М + N. (8)

а=N+1 £«) a=N+1 £«)

Здесь М и N - количество элементов на каждой из групп поверхностей. Соотношения (8) представляют собой систему линейных алгебраических уравнений относительно и (а) и . Интегралы |и*(х0,хи (х0, х^ вычисляются аналитически и образуют компоненты мат-

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

Интегралы от компонентов функций ui]■ и / вдоль граничных элементов, являющихся отрезками прямой, были вычислены аналитически. Начало отрезка обозначено точкой a(a1, a2), конец отрезка - точкой Ь(Ь1,Ь2 ). Отрезок ориентирован таким образом, что при движении от a к Ь исследуемая область находится слева, то есть граница области обходится против часовой стрелки. Результаты вычислений выглядят следующим образом:

I и11 (^, х)а5 (х)-61 C2 (d1Q1 + d 2 (Ь162 alQ3 )+ d 3Q4 d) + + d 5Q4 ^

;2 Л

d

I Ul*2 (,х ^ (x) - с | - 2d 4Й

Iи*?(х,х)(х)-с с?(а^ + d2(*102 -а^)+ d364 -d)-d461 -d564

a V

Ь

I /11 х )dS (х)-63 (б4 (04 + 1)-d 6 ),

а

| /12 ( Х )dS (х) - 63 ("6461 + d7 '1 ,

а ^ '

| Л*1 Х)dS(х) - 63 (^- + d7 '1 ,

а ^ '

Ь

I /11 Х )dS (х) - 63 (б4 (04 + 1)+ d6 ) .

2 Л

Здесь

11 8 рт(1 -п) , ^ = 3 - 4п; С3 =-4 р(1 -V),

с4 -1 - 2п

-Х25251 -Х152 + 52 57

2d51

? и2

d (51 -52 )(52 §3-§1 §4)

____ d - 5 253 -515 4 d - 515 2 (5 2 53 -5154 )

251 ’ 3 d ’ 4 d3 :

d3

d 6 =

А

25 535

34

55 5б

532 +54 51 +56

d7 - -

56

51 +56 532 +54

, d-^ +52;

61 - 1п(51 (552 +56))- 1п(5? (532 +54)), 02 - 1п(552 +56), 63 - 1п(532 +51),

04 - агС^

5153 + 5 2 54

- агСд

5155 + 5 25 6

у5154-5?53 0 ^5154-5?53 0

51 - Ь1 - а1 , 52 - Ь2 - а2 , 53 - а1 - Х1 , 54 - а2 - Х2 , 55 - Ь1 - Х1, 56 -Ь2 -Х2 , 57 - а2 Ь1 -а1 Ь2 . (9)

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

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

После решения системы и получения значений и ]а1 и /]а, перемещение в любой точке

внутренней области вычисляется по формуле (6). Деформации определяются по формуле

(х) - I \wijk (х х)А (х) - ё*(х хК(х)]*,

(10)

полученной из (6) путем дифференцирования, в силу (2). Здесь

™г]к Iх х

+а<

&

V

ёг*к(x, х)-?

Интегралы по произвольному отрезку прямой от компонент w*k и ёф также вычислены аналитически, причем результаты интегрирования выражаются через те же величины (9).

Определив деформации е ] (х) из (10), напряжения ог]- (х) можно вычислить по формуле (3).

ь

а

Б

и,/

* Л

Формулы (6), (10) и (3) справедливы для любой внутренней точки области, причем вычисление перемещений, деформаций и напряжений в одной точке не зависит от вычисления этих же величин в других точках, поэтому уровень распараллеливания для вычисления напряженно-деформированного состояния зависит только от количества процессоров.

2. Решение задач теплопроводности и диффузии. Еще одной областью применения

предложенного алгоритма являются задачи теплопроводности и диффузии, которые описыва-

ются одинаковыми уравнениями. Рассмотрим вновь некоторую двумерную область, граница которой разбита на элементы (рис. 1). Задача состоит в определении температуры (концентрации) $, которая удовлетворяет в области уравнению (для простоты изложения будем считать, что область не содержит источников)

кД$ = 0 (11)

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

на поверхности ^ : 3 = 3 *;

на поверхности Б2: д = -кпі = д*. (12)

дхг

д2 $ д2$

Здесь к - коэффициент теплопроводности (диффузии), Д$ = —— +------------2— оператор Лапласа.

дх12 дх2

Согласно методу граничных элементов [1,2], решение задачи (11), (12) можно выразить через заданные на границе значения температуры $(х ) и потока д(х), аналогично (6):

$(Х) = | [р (X, х )$* (х) - G(X, х )д(х )]& + | [р (X, х )$(х) - G(X, х )д * (х )]і5. (13)

Функции G(X, х) и р (X , х) имеют вид:

Х)=- 2-7 1п (г )> Р(Х Х) = (Х ~Хг К •

2—к 2—г

Считая вновь значения температуры и потока постоянными на каждом граничном элементе и помещенными в точку середины элемента, запишем уравнение (13) для середины каждого элемента:

N / \ N / \

$(Р) = £$(“> |р(х(р), x)ds -£ д(а) |G(х(р),х)& +

а=1 а=1 ^(а)

N + ~^М / \ N / \

+ £$(а) |р(х(р),хр-£д(а> |G(x(p),хр, р = 1..М + N• (14)

а=N+1 52а) а=1 52а)

Здесь Ми N - количество элементов на каждой из групп поверхностей. Соотношения (14) представляют собой систему линейных алгебраических уравнений относительно $(а) и д ^а).

Интегралы | р (X, х (х ) и | G(X, х ^5 (х ) суть слагаемые в уже вычисленных аналитически

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

|G(x, Х^5(х) = - + d2 (^б2-^3)+ &4 - d , Г р(X, Х^(х) = ^ .

2—к 2—

а а

3. Анализ эффективности счета задачи при замене численного интегрирования функций Грина аналитическим интегрированием. Разработанный алгоритм решения упругой и тепловой задач позволяет благодаря большому объему предварительных аналитических вычислений сократить время решения задач (а значит, увеличить число граничных элементов, на которые разбивается граница), а также повысить точность расчетов. Алгоритм позволяет распараллеливать задачу определения перемещений и тензоров напряжений и деформаций (температуры и потока тепла) за счет использования процедуры распараллеливания на трех уровнях решения данной проблемы. Во-первых, элементы матрицы системы линейных уравнений вычисляются независимо друг от друга, что позволяет расширить количество используемых процессоров до количества элементов матрицы. Во-вторых, решение системы алгебраических уравнений и нахождение искомых поверхностных напряжений и перемещений осуществляется методом Гаусса - Жордана, допускающим распараллеливание. В-третьих, компоненты перемещений и тензоров напряжений и деформаций (или значения температуры) для произвольной точки внутренней области также рассчитываются независимо, что позволяет проводить распарал-58

леливание задачи на достаточно большое количество процессоров, без больших затрат на обмен информацией между ними [3, 4]. Для проверки эффективности перехода от численного интегрирования к аналитическому был проведен расчет для задач деформирования различных двумерных областей. Собственно результаты решения никакой смысловой нагрузки не несут, поэтому здесь приводится лишь анализ скорости счета. Расчеты производились на вычислительных комплексах МВС-100 и МВС-1000.

Результаты расчетов показали, что если при приближенном интегрировании возникала ошибка в вычислении, особенно при расчете тензоров деформаций (до 10%), которую не удавалось понизить уменьшением шага интегрирования, то при счете по аналитическим формулам ошибка в вычислении перемещений, тензоров деформаций и напряжений вообще исчезала. Кроме того, получен колоссальный выигрыш во времени при вычислении перемещений и напряжений на границе области - до 3 тысяч раз, а при вычислении перемещений, деформаций и напряжений внутри области - в среднем до 600 раз.

Время решения упругой задачи складывается из 3 частей: t1 - время вычисления интегралов от компонентов функций Грина для граничных точек, 12 - время решения системы линейных алгебраических уравнений порядка М х к, где к - число точек разбиения на каждой из М частей границы области, 13 - время вычисления интегралов от функций Грина и их производных для внутренних точек области на сетке КК х КК . При замене приближенного интегрирования функций Грина на аналитическое значения t1 и t3 значительно сокращаются, что отражено в таблице. Здесь п/ р - отношение порядка системы к числу используемых процессоров, ^ , t3 - время численного интегрирования функций Грина на этапах 1 и 3, t1a , tа - время

вычисления интегралов от функций Грина по аналитическим формулам на этапах 1 и 3, к^ , к3 - соответствующие коэффициенты ускорения.

Ускорение счета при аналитическом интегрировании

п!р k, KK t[, с УЙ t J , с V t\, с Й t3 , с кУ

i 200; 4 4737 1.472 3218 418 0.220 1900

2 200; 4 2370 0.736 3220 224 0.218 1027

4 200; 4 1210 0.367 3297 108 0.246 440

8 200; 4 614 0.221 2778 59 0.328 180

16 200; 4 297 0.117 2538 29 0.393 74

32 200; 4 148 0.112 1325 16 2.047 7.8

2 500; 10 15157 5.0 3031 3292 0.815 4039

4 500; 10 7491 2.4 3121 1657 0.739 2242

8 500; 10 3746 1.2 3121 829 0.708 1171

16 500; 10 1874 0.6 3123 433 0.844 513

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

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

kDJ = 0, 0 < ф < 2р , 0 < r < R;

J = 1 + sin ф + -2 sin 3ф + cos 4ф, 0 <ф< 2р , r = R .

Аналитическое решение такой задачи выглядит следующим образом:

Г 3

J = 1 + r sin ф +—sin3ф + r4 cos4ф.

Для приближенного решения граница области - окружность с центром в начале координат и радиусом R аппроксимировалась ломаной из n отрезков прямой одинаковой длины, вершины которой лежат на окружности.

Расчеты были произведены при Я = 1. На рис. 2 показаны значения температуры внутри области, соответствующие численным решениям при п = 8 (а), п = 16 (б), п = 32 (в), а также аналитическому решению (г). На рис. 3 сравниваются значения температуры для этих четырех решений при г = 0.5 (а), г = 0.9 (б), г = 0.99 (в). По рисункам видно, что с увеличением числа граничных элементов численное решение достаточно быстро сходится к точному. При этом нет необходимости в очень мелком разбиении. График решения при п = 32 практически сливается с графиком точного решения, немного отличаясь лишь вблизи границы. Даже при п = 16 отклонение приближенного решения от точного существенно лишь при приближении к границе области, что естественно при гранично-элементном представлении.

16

3

<2 В

[и

О

-1

3

о?

и

■0 5

1 I

б

^ .6 у

3

с

Х

р# -м у

I -I

1 -1

в г

Р и с. 2. Температура внутри круговой области для разных разбиений

/1 \

"2,3,4

Чу^/

Л

Г)

/

,^-1. 5!

3,4

1,1 Ч.'1

н

ф

Р и с .3. Сравнение численных и аналитического решений тепловой задачи при различных г: а - г = 0.5 ; б - г = 0.9 ; в - г = 0.99 ; 1 - п = 8, 2 - п = 16, 3 - п = 32, 4 - аналитическое решение

а

ф

а

в

Рассмотрим теперь задачу о плоском напряженном состоянии кругового кольца при осесимметричном нагружении (рис. 4). Граничные условия имеют следующий вид:

с„

= -i^ r = r ; с r

= - P2, r = r2

22

Аналитическое решение для радиального перемещения выглядит так:

(р2 -р г г2 1 Л чрг2 ° “2

-V)

г

1

ur = — r E

-(1 + v)ip-Pr^-1 + (1 -v)

2 2 r2 - r1

Pr 2-P r

-v)_______P^ r

'1

2 2 r2 - r1

1

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

Расчеты были произведены при следующих значениях параметров: E = 210000, v = 0.3, r = 1, r2 = 2,

Pj = 100, P2 = 200. На рис. 5 показаны значения компонентов вектора перемещения внутри области, соответствующие численным решениям при n = 10 (а), n = 20 (б) и аналитическому решению (в). Из рис.5 видно, что решение внутри области даже при малом числе граничных элементов получается гладким, и лишь при приближении к границе возникают существенные погрешности. На рис. 6 показано сравнение компонентов вектора перемещения для трех полученных решений при r = 1.5 . Сравнение показывает отсутствие качественных различий численного и аналитического решений, а относительное количественное расхождение является следствием крупного разбиения.

Заключение. Разработанный параллельный алгоритм решения задач на основе метода граничных элементов позволяет быстро и с высокой точностью решать задачи с большим числом граничных элементов. В то же время, приведенные результаты решения тестовых задач дают основания полагать, что подчас нет необходимости слишком мелкого разбиения. В этом случае расчет на высокопроизводительном многопроцессорном комплексе типа МВС-1000 может быть произведен в режиме реального времени, что позволит использовать предложенный алгоритм для мониторинга прочностных свойств конструкций.

Р и с. 4. Осесимметричное нагружение кольца

ККИ

Ux П

0М15'

а начало»:

иу у -ОССй -Ofcn--ало15;

о

■now-■г,ад-

ф

CLHJ15

осот

□ □OK

Uy

у

■Dual -0 Enf.jL, D

ф

Р и с. 5. Перемещения внутри кольца при осесимметричном нагружении

а

и

x

б

Р и с. 6. Сравнение численных и аналитического решений задачи о нагружении кольца:

1 - n = 10, 2 - n = 20, 3 - аналитическое решение

БИБЛИОГРАФИЧЕСКИЙ СПИСОК

1. Brebbia C.A., Telles J.C.F., WrobelL.C. Boundary Element Techniques. Springer-Verlag Berlin Neidelberg New-York Tokyo, 1984. 524 p.

2. Бенерджи П., Баттерфилд Р. Методы граничных элементов. М.: Мир, 1984, 494 с.

3. Гасилов В.Л., Думшева Е.С., Зенкова Е.С., Федотов В.П. Численное моделирование упругой задачи на многопроцессорных вычислительных системах // Алгоритмы и программные средства параллельных вычислений. Екатеринбург: УрО РАН. Вып. 6. 2002. С. 104-124.

4. Думшева Е.С., Зенкова Е.С., Федотов В.П., Спевак Л.Ф., Привалова В.В. Численно-аналитический алгоритм для решения задач упругости, теплопроводности, диффузии // Алгоритмы и программные средства параллельных вычислений. Екатеринбург: УрО РАН. Вып. 7. 2003. С. 70-86.

Работа выполнена при поддержке гранта РФФИ №04-01-00274.

Поступила 28.05.2004 г.