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



УДК 519.632.4

A.B. Пашковский

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

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

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

Основная идея метода стандартных элементов (МСЭ) состоит в заполнении расчётной области совокупностью стандартных элементов (СЭ), в которых известны аналитические решения соответствующих краевых задач, найденные классическими методами. Такой подход позволяет не только резко сократить степень дискретизации расчётной области, но и на порядки повысить точность решения. Метод стандартных элементов, впервые введенный в [1] для тепловых расчётов электрических машин, при заполнении расчётных областей предполагал исполь-

зование таких СЭ, как: круг, сектор, сектор кольца и др. Комбинированное использование МСЭ с МКЭ также отличалось резким повышением точности расчётов.

Дальнейший анализ геометрии электромагнитных устройств (ЭУ), а также разнообразных КОС, показал значительные перспективы использования в моделировании электромагнитных и температурных полей СЭ типа: прямоугольник и параллелепипед. Результаты, полученные в [2, 3], показали стройность лежащей в основе метода математической теории и соответствующих выводов. Данные [4], кроме того, проиллюстрировали достаточную простоту использования СЭ прямоугольной геометрии для заполнения двух-и трехмерных расчётных областей, удобную склейку СЭ межцу собой, возможность комбинирования подобных СЭ с конечными элементами в расчётной области при сохранении высокой точности расчётов.

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

Первоначально для сравнительной оценки МСЭ и МКЭ получены распределения магнитного потенциала в однородной области (рис. 1). Геометрия области аналогична замкнутой системе П-образного магнита. Зазор в магнитопро-воде принимался равным 2 мм, что соответствовало размерам в реальной системе. Ставилась

^Научно-технические ведомости СПбГПУ 6'

задача не только получить численные решения краевой задачи МСЭ и МКЭ, но и сопоставить полученные данные с её аналитическим решением. В качестве аналитического решения в расчётной области в порядке усложнения рассматривались модельные функции (рис. 2)

щ = 7 + Зз> + 4х + 5хз'; Щ = У3 - Зх2)';

щ = 51П(у _ 9 )8т(х)/10 + 1, (1)

определяющие граничные условия (2) в задаче Дирихле для уравнения Лапласа:

сКу(1/мвга£1£/) = 0,

Мс, = иь ЦС2 = иъ Цсз = и3, = Щ. (2)

Также для существования решения в аналитическом виде, считалось, что магнитные проницаемости материалов области одинаковы и

источники поля в области отсутствуют. Для ис-▲

У

4Х 47 | 4Г, | 4^ | 44 | 43

17 39 40 4] 4?

Зь 35 34 • 33 32 31

.и }(■ 71 ?8 3"

24 23 22 21 20 19

П 14 И V. .17 ' ' 1»

12 11 10 • 9 7

1 2 3 4 5 Л

-►

Рис. 1. Расчётная область (48 СЭ)

пользования МСЭ расчётная область была разбита на 48 СЭ (см. рис. 1) в виде прямоугольников, что продиктовано удобством их формы для её заполнения. Для контроля точности было выделено девять точек (см. рис. 1), включая точки в угловой части центрального сердечника и в узком зазоре, рассчитываемые обычно с наименьшей точностью. Значения аналитического решения в точках находились непосредственным вычислением значений тестовой функции. Для получения картины двухмерного поля МКЭ применялись программные пакеты БЕММ и МагЫаЬ. Разбиение расчётной области на конечные элементы представлено нарис. 3.

Полученные результаты показали, что МСЭ позволяет резко сократить степень дискретизации расчётной области (табл. 1).

Очевидны преимущества МСЭ и в точности расчёта поля при наличии в области узких включений и быстро осциллирующих функций.

Существует мнение, что увеличение степени дискретизации расчётной области при использовании МКЭ всегда влечет повышение точности. Как показали численные эксперименты,

Рис. 3. Разбиение расчётной области на КЭ (3784)

Таблица 1 Результаты расчёта значений тестовых функций

Метод Число Число IV, »2 Из

решения узлов элементов ^тах? /о ^гпах? ^шах?

Аналитика 0 0 0 0 0

МКЭ 1954 3784 0,0560 1,3 4,21

МСЭ 63 48 0,0612 0,1330 0,6035

проведенные нарассмотренном примере, увеличение числа элементов с 3760 до 15136, а затем и до 60544 (более чем в 16 раз) не принесло значительного уменьшения величины погрешности МКЭ (табл. 2).

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

Однако наибольший интерес представляет расчёт областей с кусочно-однородной структурой. Рассмотрим результаты численных экспериментов (табл. 3) по решению краевой задачи (2) в расчётной области (см. рис. 1) с учётом реальных магнитных характеристик входящих в нее материалов (цвоздуха = 1, цмеди = 1, цсталим-45 = 4985) и тестовых функций вида (1).

В силу отсутствия аналитического решения задачи Дирихле для уравнения Лапласа при кусочно-однородной структуре области возьмём за эталон значения, полученные по МКЭ. Следует заметить, что в силу полученных выше фактов такое предположение является спорным. Тем не менее становится очевидным факт увеличения разницы в значениях сравниваемых методов при

Таблица 2

Результаты расчёта МКЭ значений И^ при увеличении степени дискретизации

3760 элементов ^тах, 0 1,3 1,2 1,2 1,2 1,1 1,24 1,1 1,2

15136 элементов ^тах, 0 1,39 1,24 1,23 1,32 1,31 1,29 1,24 1,25

60544 элементов ^тах, 0 1,44 1,26 1,26 1,2 1,28 1,31 1,34 1,33

Таблица 3 Результаты расчёта для тестовых функций

Метод решения Число узлов Число элементов И7! ^тах, 0 ИА ^тах, 0 ^тах, 0

МКЭ 1954 3784 - - -

МСЭ 63 48 0,064 2,6 20

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

На основании вышеизложенных результатов численных экспериментов, их сравнения с аналитическими решениями, результатами, полученными в [1], можно сделать следующие выводы:

МСЭФ позволяет резко сократить объём работы по дискретизации области;

МСЭФ обеспечивает более высокую точность расчёта однородных областей с осциллирующими функциями и наличием узких включений, чем МКЭ;

необходимо разработать математическую модель кусочно-однородной среды для полноценной сравнительной оценки точности МСЭ и МКЭ.

Математическая модель КОС

Рассмотрим несколько прямоугольных СЭ (рис. 4).

Зададим в каждом прямоугольнике гармонические модельные функции соответственно их, и2,... вида:

их = к/хж\(ких-у)й\(ких-х), и2 = к/2ж\(ких-у)с\\(ких-х)... ин = к/$со$,(ких-у)й\(ких-х), и1 = к/7со?,(ки{-у)с\\(ки{-х)...

Пусть в области поставлена краевая задача Дирихле для уравнения Лапласа:

сЦ¥(1/цёгаа£0 = 0, 2

где условия^ вводятся в соответствии с выбранными в СЭ модельными функциями.

Проведём первоначальные рассуждения для пары функций их, и2. Функция решения в области должна бытьнепрерывна, а следовательно, на

ы

ь\

8 7 6 5

1 2 3 4

а 1

а.2

дЗ

а4

Рис. 4. Кусочно-однородная область

Научно-технические ведомости СПбГПУ 6'2009

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

и 11л- = й2 = и2\х = а2- мз1.г = йЗ = и2\х = йЗ О)

И Т. Д.

Таким образом, связь между коэффициентами, например на границе я2, с использованием (3) определена:

к/2 = к/х-1к{а2-ких). (4)

Рассмотрим поведение производных функций на этой же границе:

1Д/2М1(ких ■ а2)и2Хга2 = = {/к/х/сЪ(ких ■ а2)и1х\х=а2. (5)

С учётом стандартного вида условий на внутренних границах расчётной области в электромагнитных расчётах и (5), можно считать, что в первом СЭ задана среда с щ = к/1- сЬ(кщ ■ а2), во втором — с ц2 = к/2 ■ яЬ(ки{ ■ а2) и ц2 = М-1 • Ж2(ки1 ■ а2).

Рассуждая аналогично, рассмотрим "склейку" 1-го и 8-го СЭ вдоль оси оу:

к/х ■ &к(ких ■ х^пт^А:^ • Ь2) = = к/ц ■ &к(ких ■ х)со&(ких ■ Ь2)\ -1/к/н/ш1(ких ■ Ь2)иЬу\х=Ь2 = = 1/к/х/со$(ких ■ Ь2)и1у\у=Ь2.

Ранее было получено, что в первом прямоугольном СЭ среда с ц, = к/х ■ ch(^«l • а2), тогда

-\/к^т{ких ■ Ь2)/сЩих ■ а2)сс&(ких ■ Ь2)иНу\у = Ь2 = = 1/к/х/сЬ(ких ■ а2)и1у\у = ь2-ц8 = —• tg(A;м1 • Ь2) ■ сЪ.(ких ■ а2) и ц = -щ-Щ2{ких ■ Ь2).

Таким образом, алгоритм определения связей коэффициентов и магнитных проницаемос-тей СЭ как по горизонтали, так и по вертикали подобны. Обобщая данную методику на рассмотренную выше область из сорока восьми СЭ и используя формулы типа (4), (5), получаем шахматную структуру неоднородности расчётной области (рис. 5). Таким образом, построена модельная краевая задачадля кусочно-однород-

к

У

0 к

Рис. 5. Структура неоднородности

ной среды с известным в каждом СЭ аналитическим решением.

Погрешности в расчёте потенциала для девяти контрольных точек, включая точки зазора, приведены в табл. 4.

Проведенные исследования показали, что МСЭ позволяет:

резко сократить объём работы по дискретизации области при сохранении точности расчётов благодаря выделению в расчётной области стандартных элементов, аналитическое решение поставленной задачи в которых известно;

достичь высокой точности расчётов благодаря замене аппроксимации расчётной функции полиномами (как в МКЭ) на аналитическое представление функции и получать непрерывно дифференцируемые решения внутри стандартных областей, что важно для вычисления потоков;

использовать классические решения для учёта особенностей поведения решения в окрестности угловых точек области;

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

Таблица 4

Значения погрешностей методов

№ п/п 1 2 4 5 6 7 8 9

МКЭ 3,159 6,503 0,0071 0,0137 0,2380 0,0067 6,259 3,004 2,735

МСЭ (> 4раз) 0,6092 1,627 0,0058 0,0063 0,0920 0,0062 1,546 0,7485 0,6763

Дополнительно к выводам, сделанным выше, можно добавить:

МСЭ обеспечивает достаточную стабильность точности расчёта при усложнении вида функции и росте значений её градиента в расчётной области.

Рассмотренные модельные функции изменяются достаточно плавно, что приводит к стабильной точности как МКЭ, так и МСЭ.

Преимущества МСЭ по точности проявляются для быстро изменяющихся или быстро осциллирующих функций.

СПИСОК ЛИТЕРАТУРЫ

1. Пашковский A.B. Расчет температурных полей в элементах электрических машин и аппаратов: Дис. ... канд. техн. наук / НГТУ (НПИ). Новочеркасск. 1989. 216 с.

2. Пашковский A.B., Пашковская И.В. Прямоугольный стандартный элемент в моделировании температурных и электромагнитных полей в кусочно-однородных средах // Электромеханика. 2003. № 3. С. 9—12.

3. Пашковский A.B., Пашковская И.В. Прямоугольный стандартный элемент с условиями Ди-

№

С. 71-73.

4. Пашковский A.B., Пашковская И.В. "Склеенные" прямоугольные стандартные элементы в решении полевой задачи // Электромеханика.

№

УДК 622.24

A.A. Цуприков, В.Г. Чередниченко

РАСЧЕТ ДИНАМИЧЕСКИХ НАГРУЗОК В БУРИЛЬНОИ КОЛОННЕ

ПРИ ПОСАДКЕ НА ЗАБОЙ

Однойизхаракгфисппссистемы управления технологическим процессом (ТП) роторного бурения нефтяных и газовых скважин считается частота опроса первичных датчиков параметров контроля иуп-равления. Интервал между считываниями данных во многом определяется временем переходных процессов (ПП), происходящих в скважине в результате изменения режимных параметров. Основной параметр управления для ТП бурения — осевая нагрузка надолото, которая на60—70 % обеспечивает эффективность углубления скважины и является главной причиной возникновения ПП в скважине.

Вывод математической модели ПП, возникающих в скважине в результате изменения нагрузки на долото, приведен в настоящем сборнике (см. статью В.Г. Чередниченко, A.A. Цуприкова "Моделирование динамических нагрузок в бурильной колонне при посадке на забой"). В указанной статье приводятся результаты моделирования ПП в среде программирования Matlab 7.5 (R2007b)*.

* Алексеев Е.Р., Чеснокова O.B. MATLAB 7. М.: НТ Пресс. 2006. 464 е.: ил. (Самоучитель.)

Результаты вычислений

Программные вычисления по модели ПП проводились в безразмерном виде при следующих значениях параметров: ук=1,5;а=1,5;?=10; к = 200, где ук — скорость движения колонны

при посадке; а = —

элементарного сечения 5 колонны труб с плотностью материала труб р; р — коэффициент, пропорциональный вязкости бурового раствора с р

ни; к — количество сечений колонны, в которых просчитывается деформация.

Профили деформации колонны по её длине (в относительных единицах деформации) враз-личные моменты в течение 10 относительных единиц времени представлены на рис. 1. Задняя плоскость соответствует нижней части бурильной колонны, т. е. забою (х= 0), передняя — верхней части, т. е. устью (х = 1).

На рис. 2 показано изменение во времени относительной деформации элементарного участка колонны сечением 5 и длиной ¿/хна концах