CC BY
7
УДК 621.372.001
ВЕРИФИКАЦИЯ СХЕМОТЕХНИЧЕСКИХ МОДЕЛЕЙ В МЕТОДЕ КОНЕЧНЫХ ЭЛЕМЕНТОВ
В.Ф. Булавин
Проанализирована концепция континуальности физико-математической модели КЭ, которое требует сравнительной оценки на контрольном примере схемотехнической и классической реализаций МКЭ и анализа итоговых результатов. В качестве контрольного материала использована тестовая задача, имеющая аналитическое решение. Для этого примера приведены результаты, иллюстрирующие классическую технологию МКЭ, и реализованы альтернативные подходы. Вычислительный эксперимент преследует цель всесторонней проверки альтернативных моделей КЭ для всех типов шаблонов. Анализ результатов по схемотехнической и классической реализациям проводится по двум критериям: близость к аналитическим значениям в контрольном примере и обусловленность глобальной матрицы жесткости. Затронут вопрос рационального выбора пробной функции при расчете градиентов поля.
Ключевые слова: формула Грина, метод конечных элементов, схемотехническая модель, матрица жесткости, аппроксимация, обусловленность алгебраических уравнений.
Классическая парадигма определяет математическую модель (ММ) конечного элемента (КЭ) в виде матрицы жесткости (МЖ) как единственную [1 - 4]. Схемотехническая концепция [5 - 7], в основе которой лежит метод интегральных тождеств [8], оперируя понятием континуальности, позволяет расширить круг возможностей. Консервативность модели обеспечивается использованием интегроинтерполяционным методом на базе инвариантных базисных функций [9].
Алгоритм построения контрольных сечений и теоретическое обоснование при формировании семейства схемотехнической модели КЭ изложен в [5].
Вызовы о континуальности модели КЭ требуют сравнительной оценки на тестовом примере схемотехнической и классической реализаций МКЭ и анализа итоговых результатов. В качестве контрольного материала использована классическая задача Сен-Венана [10] о кручении цилиндрического стержня некругового сечения, имеющая аналитическое решение. В [4] для этого примера приведена "ручная" технология, иллюстрирующая классическую реализацию FEM.
Расчетная модель в виде 1/8 части исходной области и варианты конечно-элементной сетки (виде сплошной и пунктирной линий) показаны на рис. 1, а, б.
Нестандартное разбиение преследует цель всестороннего вычислительного эксперимента, проверки континуальности схемотехнической модели КЭ и альтернативных возможностей для всех типов шаблонов. В каждом варианте элементные МЖ определялись для всех точек, соответствующих предложенным схемотехническим решениям.
В расчетной области поля функция напряжений ф (х,у) подчиняется дифференциальному уравнению
Э 2ф Э 2ф
+ 200 = 0; ф Эд = 0,
Эх2 Эу2
Р 1
где G=0,8x107 [H/Cм2]- модуль сдвига материала; ^щ) [рад/См] -
угол закручивания на единицу длины. Размеры образца и исходные данные взяты из [4].
Результаты расчетов на стандартной конечно-элементной сетке (сплошная линия) приведены в табл.1.
Рис. 1. Расчетная модель контрольного примера и два варианта
разбиения на конечные элементы
Комментарии к результатам табл.1. Различие результатов в столбцах 2 и 3 носит принципиальный характер. При одинаковых элементных и глобальной МЖ имеет место несовпадение правых частей конечно-элементной системы уравнений. Причина этого кроется в принципе их формирования. Для каждого КЭ автор [4] распределяет источники поля по
узлам как 2G 9 ^ = 29,07, т. е. без учета реального положения точки О -
центра описанной окружности, относительно которой формируется МЖ элементов.
На рис.2 на примере первого фрагмента показано, что распределение источников поля по узлам (или в схемотехнической интерпретации определение эквивалентных источников тока в соответствующей аналоговой модели), с учетом положения точки О, осуществляется пропорционально площадям сечений Ае2 и Ае4, а именно: для узла 1 -
е $
Зх = 20$-^ = 21,817; для узла 2 - /2 = 200-^ = 43,633 и узла 4 -
е
Л = 200-^ = 21,817. 4 4
Таблица 1
Результаты сравнительного_расчета тестовой задачи (вариант 1)
Номер узловой точки Точное решение Приближенное решение [4] Решение при построении МЖ относительно точки О-центра описанной окружности Решение при построении МЖ относительно точки Р- центра окружности Эйлера Решение при построении МЖ относительно точки О-центра тяжести Решение при построении МЖ относительно точки Н-ортоцентра Решение при построении МЖ относительно точки Q-центра вписанной окружности
1 2 3 4 5 6 7
Ф1 205,6 218,04 196,4 218,12 210,12 - 213,62
ф2 160,0 159,9 152,72 141,8 147,27 - 144,9
ф4 126,4 123,56 119,99 141,8 132,32 - 136,44
Число обусловленности глобальной МЖ (матрицы узловых п] роводимостей)
- 10,89 10,89 6,31 7,51 - 6,9
Каковы бы ни были доводы в пользу классического подхода и как бы убедительно они ни звучали, их нельзя противопоставить интегро-интерполяционному методу и фундаментальному закону-постулату Максвелла, которые требуют иного решения, приводя, к консервативному результату.
Результаты вычислений в 6 столбце отсутствуют в виду некорректности соответствующей схемотехнической модели (ортоцентр, точка Н, всех КЭ совпадает с узловой точкой).
На рис. 1, б пунктиром представлено альтернативное разбиение, особенность которого заключается в наличии трех тупоугольных фрагментов (элементов 1,3 и 4). Вызовы статьи не позволяют использовать классическую МЖ для указанных КЭ, а, следовательно, провести расчет в целом. И хотя данная ситуация (с такой долей "плохих" элементов) носит искусственный характер, автор тем самым стремится показать, что наличие тупоугольных фрагментов не является препятствием для проведения расчетов, при условии альтернативных решений для МЖ.
А Л
л / уУ
О ¿Г^ Л
Л®
/ (¿) 5: \
1 . / -Д
Рис. 2. Схемотехническая модель прямоугольного шаблона с центральной точкой контрольных сечений в центре описанной
окружности
359
Результаты вычислительного эксперимента на нестандартной сетке сведены в табл. 2.
Таблица 2
Результаты сравнительного расчета тестовой задачи (вариант 2)
к и сг1 о
н «
о
(Я
о ч
СО
а <и
5! о
к
<и к
X <и
а
<и а <и о X сг1 о н
1
<и к
X <и
а
<и а <и
О ,—] X ^ И ^
<и *
X
ч ю к
£
2
к я
5 О
о
к
и ^
о н о
д ч <и <и
а н о о с
к
X
<и 0
<и Рч
к н о о X
а &
о «
о X X а о к с о
<и к
X <и
д
<и Рч
р ч
а 5®
5 е
о
о X
>у &
О
оЗ р
<и я
4
I—г
к я
5 6
° а
^ В
о И
о 2
X Н
к 8
<и <и
К Н
И К
5
к к
я ,
о
р
н о о с
<и К X <и
д
<и Рч
ОЗ р
Н
и а
Я о о
к к
° а ^
а й ^ нал У К *
В £ 1:5 « § °
р ч 5
<и <и
К Н Г1 К о о X н о
X <и
д
<и Рч
О
К Сч (Я
7
ф
205,6
132,02
253,6
184,77
197,04
ф2
176,9
92,16
192,89
137,37
147,1
ф4
126,4
65,68
187,04
120,64
137,94
Число обусловленности глобальной МЖ (матрицы узловых проводимостей)
2,09
6,18
5,74
5,85
3
6
Комментарии к результатам табл.2. Решения в столбце 3 отсутствует ввиду некорректности элементных МЖ для фрагментов 1, 3 и 4 (точка О -вне контура КЭ), а также в столбце 6-ввиду некорректности элементных МЖ для всех шаблонов.
Результат в столбце 2 получен по классической технологии МКЭ при формализации процесса вычислений (не обращая внимания на отсутствие физического смысла у элементных МЖ фрагментов 1, 3 и 4). Приведем глобальную МЖ в этом случае (с учетом граничных условий):
' 0,65 - 0,75 0,1 Л
К= - 0,75 2,292 - 0,292 ч 0,1 - 0,292 1,342
V ' ' ' У
Наличие положительных недиагональных элементов является причиной того, что такой структуре невозможно сопоставить физический аналог в виде схемотехнической модели, реализованной на пассивных элементах. Этот результат вступает в противоречие, как на этапах формирования, так и окончательном виде с формулой Грина и постулатом Максвелла [11-13].
Методология распределения источников поля по узлам в классической модели для рассматриваемого случая осуществляется как
20 = 29,07,
что после ансамблирования и учета граничных условий
дает результат для правой части СЛАУ в виде: В=[23,26; 93,04; 74,43]т. Этот ответ вступает в серьёзное противоречие с альтернативными решениями, а по сути является некорректным.
360
Оставаясь на позициях трактовки по Грину схемоанализа и МКЭ [5-7,14,15], ясно, что такое условие выполняется только для равностороннего КЭ. По мере искажения формы шаблона от указанной наблюдается рост методологической погрешности в итоговом результате. Указанные факторы объясняют слабый (некорректный) итог во 2-м столбце.
Помимо требований близости численного решения к аналитическому (столбец 1), в обоих вариантах в качестве критерия для сравнении вычислительных свойств глобальных МЖ (матриц узловых проводимостей)
1
приводятся их числа обусловленности ( max ) (1- собственное число МЖ).
1 min
Из табл.1 и табл.2 видно, что альтернативные решения обладают заметно лучшими показателями.
В табл.2 аналогичные данные во 2 столбце не обсуждается, как полученные для матрицы, не имеющей физического смысла, и приводится здесь по принципу "as is". Но в контрольных решениях обусловленность даже лучше, чем в аналогичных данных из табл.1
В целом, как указывает автор [4, стр.92], "четырех элементов мало для получения приемлемой точности решения, но достаточно для иллюстрации".
Целесообразность применения комбинированного алгоритма, когда для "хороших" КЭ используется классическая МЖ, а для "тонких" фрагментов один из альтернативных вариантов возможна, но не очевидна. Такая стратегия потребует включения дополнительного логического блока, что нарушит однородность вычислительной схемы.
Ансамблирование матрицы узловых потенциалов (глобальной МЖ) при схемотехническом подходе может быть осуществлено методами теории электрических цепей. Методологически эта процедура выглядит как применение первого закона Кирхгофа для всех внутренних узлов расчетной области. В узловом анализе электрических цепей имеются эффективные алгоритмы формирования уравнений, например, метод поэлементного вклада, позволяющий сформировать требуемую матрицу проводимостей [12]. Отметим, что такой важный показатель для глобальной МЖ, как степень разреженности, не ухудшается.
В процессе составления узловых уравнений следующая ситуация должна быть принята во внимание. Если узел с номером "i" является внутренним, т.е. окружен тремя или более КЭ, то в этом случае группа слагае-
мых вида -ц
ЭфЛ
Эп
• x +
Эп
V
•(a-ha )
и аналогичных, соответ-
ствующих всем ребрам, примыкающим к рассматриваемому узлу, в сумме обращается в ноль. Эквивалентный источник тока в рассматриваемом узле образуется как сумма слагаемых вида X Рк, обеспечивая баланс потоков между КЭ. Реализация в граничных узлах краевых условий всех видов не вызывает сложностей.
с
a
Выбор пробной функции для расчета градиентов поля. Декомпозиция области поля с границами сложной геометрии приводит к образованию КЭ разнообразного вида. Для остроугольных и прямоугольных фрагментов инвариантные решения фк равноценны, но использование тупоугольных
шаблонов практика прикладных расчетов не рекомендует [2-4], объясняя это большой погрешностью. Однако именно такие шаблоны при аппроксимации криволинейных границ чаще всего имеют место.
Вопрос равноценности базисных функций важен при расчете градиентов поля. Конечно-элементная практика таких вычислений строится на базе только одной из трех возможных, составляющих линейную оболочку, функций фк. Опишем простой алгоритм, при выполнении которого эта
проблема легко решается путем рационального выбора базисной функции КЭ [14,15].
На рис. 3, в местной системе координат (х,у), ассоциированной с ребром шаблона "с", линейную пробную функцию запишем в форме
ф(х,у):
1
х
С2 X у
с
К х с у
Хфг. +
х с1 х у
Кс Х с у
Хф / +
/ Л
У_
V Кс У
Хф
т
ф т - Ф
К
* Л
У +
Ф / -ф,
• х + ф,-
где принято ф
Ф / ф
с
• с1+ф,
Однако значение ф лежит вне контура шаблона и не является элементом линейной оболочки рассматриваемого КЭ. Эта величина может быть найдена только из соседнего шаблона. Можно утверждать, что эта аппроксимация, удовлетворяя исходному уравнению и краевым условиям, не является при этом решением задачи о распределении полевой функции на КЭ. Следовательно, в системе координат (х, у) при описании состояния КЭ не рекомендуется пользоваться финитной функцией симплекс-элемента для сохранения заданного уровня погрешности [14,15]. Таким образом, указанный выбор частного решения некорректен, так как приводит к неконтролируемому результату.
Ситуация разрешима, если воспользоваться другим частным решением в локальной системе координат, связанной со стороной "а" (рис.4). В указанной ситуации в осях (л, X) базисная функция выглядит следующим образом:
г *\
ф ]-ф Ь а
Х+
ф! -фт
а
•Л + фт
ф
ф! -фт
а
• а1 + фт.
В последнем случае все величины, входящие в определение базисной функции, принадлежат линейной оболочке треугольника.
с
с
*
С2 С1
Рис. 3. Некорректный выбор системы координат и пробной функции
4 Ф/
ЪуГ
ha
Ф* N. Ф'
фш <71 «2
(7 1
Рис. 4. Корректный выбор системы координат и пробной функции
Наличие свободного заряда требует детализация составляющих поля. Так, в этом примере после коррекции окончательный результат для градиентов в направлении осей получаем в виде
дФ _ 1 f a1 m + m a2 m 1 XP • ЭФ _ ji -jm
1--ji +jj--jm \ ; .
dq ha у a a J 2ц дц a
При необходимости следует учесть точечные источники, что было рассмотрено в первой части данной работы.
Сказанное выше является одной из причин, препятствующей использованию тупоугольных фрагментов в классической интерпретации МКЭ. Таким образом, указанный выбор частного решения приводит к контролю заданного методологического уровня погрешности.
Окончательно требования к вычислению градиентов поля в тупоугольных фрагментах в МКЭ можно сформулировать следующим образом [14,15]:
для тупоугольных КЭ описание состояния элемента осуществляется базисной функцией, которая должна быть определена в локальной системе координат, связанной с самым длинным ребром:
необходимые градиенты поля следует вычислять на основе базисной функции, ассоциированной с этой стороной КЭ;
переход к числовым значениям в глобальной системе координат осуществляется на основе стандартной процедуры.
Назначенные в качестве центров формирования сечений две точки: центр тяжести КЭ и центр окружности, вписанной в шаблон, обеспечат результат с заданной методологической погрешностью. Это позволит в процессе реализации МКЭ снять все ограничения на вид КЭ и тем самым расширить его возможности. Методологическая погрешность решения при выполнении указанных требований всегда будет сохраняться на уровне
o(h2), где h _ max(a,b,c).
Выводы
1. Интегроинтерполяционный метод позволяет сформировать для симплекс-элементов семейство схемотехнических моделей, часть из которых имеет нестандартные матрицы жесткости. Число возможностей составляет множество, мощность которого континуум.
2. На контрольном примере проведена верификация концепции о континуальности схемотехнической модели конечного элемента. Оценка полученных численных решений проводится по двум критериям: близость узловых переменных к аналитическим значениям и обусловленность глобальной матрицы жесткости.
3. Результат проверки показывает, что обусловленность конечно-элементных уравнений для альтернативных решений значительно улучшается.
4. При схемотехнической реализации МКЭ учет краевых условий любого вида сводится к расчету эквивалентных источников тока в граничных узлах, что автоматически может быть принято во внимание в алгоритме формирования конечно-элементной системы уравнений на базе теории цепей.
5. Указаны универсальные решения для элементных матриц жесткости, реализующие методологическую погрешность метода конечных элементов, для произвольной по формы вычислительных шаблонов.
6. Вычисление градиентов поля в тупоугольных линейных КЭ требует целенаправленного выбора базисной функции. Использование таких фрагментов упрощает стратегию автоматического разбиения области поля на КЭ, позволяя точнее аппроксимировать границы сложной конфигурации.
Список литературы
1. Шакиров М.А. Теоретические основы электротехники. Новые идеи и принципы. Схемоанализ и диакоптика. Декомпозиционные алгоритмы. СПб.: Изд-во С.-Петербургского университета, 2001. 212 с.
2. Зенкевич О. Конечные элементы и аппроксимация. М.: Мир, 1986. 318 с.
3. Ильин В. П. Методы и технологии конечных элементов. Новосибирск: Изд-во Института вычислительной математики и математической геофизики Сибирского отделения Российской академии наук, 2007. 371 с.
4. Сегерлинд Л. Применение метода конечных элементов. М.: Мир, 1979. 392 с.
5. Булавин В.Ф. Континуальные схемотехнические модели в методе конечных элементов // Электричество. 2015. №1. С. 39-51.
6. Булавин В.Ф., Благовестова М.Е. Метод конечных элементов: альтернативные решения для матрицы жесткости // Известия вузов Электромеханика. 2017. Т. 60. №1. С. 5-13.
364
7. Булавин В.Ф., Булавина Т.Г., Яхричев В.В. Инженерный анализ и новые технологии в методе конечных элементов // Фундаментальные и прикладные проблемы техники и технологии. 2018. № 2(328). С. 109-120.
8. Смирнов В.И. Курс высшей математики. Т. II. 24-е изд. СПб.: БХВ-Петербург, 2008. 848 с.
9. Дружинин Г.В. Построение базисных функций и их применение к краевым задачам механики сплошной среды // Прикладная механика и техническая физика. 2003. Т.44. №6 (262). С. 35-43.
10. Канторович Л.В., Крылов В.И. Приближенные методы высшего анализа. М.: Изд-во физ.-мат. лит., 1962. 709 с.
11. Смирнов В.И. Курс высшей математики. СПб.: БХВ-Петербург, 2010. Т. 3. Ч. 1. 400 с.
12. Демирчян К.С., Нейман Л.Р., Коровкин Н.В и др. Теоретические основы электротехники: учебник для вузов. СПб.: Питер, 2004. 463 с.
13. Киреев В.И., Пантелеев А.В. Численные методы в примерах и задачах. М.: Высш. шк., 2008. 480 с.
14. Булавин В.Ф. Метод конечных элементов: возможность применения симплекс-элементов произвольной формы // Информатизация процессов формирования открытых систем на основе СУБД, САПР, АСНИ и систем искусственного интеллекта: материалы 7-й МНТК. Вологда: ВоГ-ТУ, 2013. 243 с.
15. Булавин В.Ф. Метод конечных элементов и выбор линейной базисной функции. Автоматизация и энергосбережение машиностроительного и металлургического производств, технология и надежность машин, приборов и оборудования // Материалы 8-й МНТК. Вологда: ВоГТУ, 2013. 264 с.
Булавин Вячеслав Федорович, канд. техн. наук, доцент, bulavin35@,mail. ru, Россия, Вологда, Вологодский государственный университет
VERIFICA TION OF SCHEMO TECHNICAL MODELS IN THE FINITE ELEMENT METHOD
V.F. Bulavin
The concept of continuity of the physical and mathematical model of FE is analyzed. It requires a comparative assessment on a test example of the circuitry and classical implementations of the FEM and analysis of the final results. The test task having the analytical solution is used as the control material. For this example, the results illustrating the classic FEM technology are given and alternative approaches are implemented. The computational experiment is aimed at comprehensive verification of alternative FEM models for all types of templates. The analysis of results on schematic and classical implementations is carried out on two criteria: closeness to analytical values in a control example and conditionality of a global matrix of rigidity. The question of rational selection of the trial function when calculating the gradients of the field will be raised.
Key words: Green formula, finite element method, schematic model, the stiffness matrix, the approximation, conditionality of algebraic equations.
Bulavin Vyacheslav Fedorovich, candidate of technical sciences, docent, bulavin35amail. ru, Russia, Vologda, Vologda State University
365
