разное

Метод конечных элементов для решения дифференциальных и интегральных уравнений

Метод конечных элементов (МКЭ, по-английски FEM — Finite Element Method) состоит в том, что все тело разбивают на некоторое коли­чество частей (элементов) конечного (а не бесконечно малого) объема, на­столько простых по форме и внутреннему устройству, что интегрирование уравнений внутри каждого из них не вызывает затруднений. Для прибли­женного решения задачи в целом необходимо обеспечить только стыковку элементов между собой. Условия стыковки записывают в виде алгебраиче­ских уравнений. В некоторых случаях эти уравнения независимы друг от друга и могут быть решены по одиночке, но чаще они образуют систему уравнений, порядок которой зависит от числа конечных элементов.

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

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

Таким образом, наличие погрешности при использовании конечного элемента не является свидетельством его непригодности, если эта погреш­ность убывает при уменьшении размеров элемента и в пределе стремится к нулю. Это не зна­чит, что все виды элементов равноценны. При прочих рав­ных условиях следует предпо­читать элементы, дающие на­именьшую погрешность при данных размерах и более бы­строе ее убывание при умень­шении размеров. Рис. 2.2. Виды конечных элементов

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

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

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

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

В мире существуют десятки коммерческих программных комплек­сов МКЭ (наиболее известны NASTRAN, ANSYS), в том числе специали­зированных для решения сварочных и других технологических задач (SYSWELD, MARC). На их создание были затрачены большие усилия, и тем не менее нельзя назвать ни одного, пригодного для решения всех возникающих задач. Поскольку внутренняя часть программного ком­плекса является «черным ящиком», то довольно трудно бывает самостоя­тельно приспособить его к решению задачи, не предусмотренной разра­ботчиками. Еще одно обстоятельство оказывается не в пользу коммерчес­ких комплексов — огромные размеры текста программ (сотни мегабай­тов) и обязанность поддерживать совместимость новых версий с преды­дущими. Теряется гибкость и с трудом осваивается решение принципи­ально новых задач, а изменения в аппаратном обеспечении (например,

переход с больших ЭВМ на персональные компьютеры) имеют для таких комплексов катастрофические последствия.

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

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

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

Конечный элемент (КЭ) позволяет установить, какой ток потечет из клетки 2 в клетку 1 через их границу s, если известна разность потенциалов между этими точками. Распределение потенциала по длине элемента ап­проксимируем полиномом первой степени, напряженность поля во всех

точках элемента одинакова и равна Е-^г^х. Плотность тока, согласно £

закону Ома (2.1), j = —, поскольку материал изотропный и сопротивление Р

во всех направлениях одинаково. Ток через элемент равен

и2-их

Таким образом, КЭ эквивалентен сопротивлению R, включенному между узлами 7 и 2. Оно равно сопротивлению проводника из такого же ма­териала, длиной / и сечением s (см. рис. 2.3). Если добавить КЭ для каждой пары соседних узлов, то они покроют всю пластину. Если считать, что все заряды, попавшие в одну из клеток, находятся в узле, то всю пластину мож­но изобразить в виде принципиальной электрической схемы (рис. 2.4). Со­гласно классификации КЭ, рассмотренный элемент следует называть стерж­невым линейным двухузловым.

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

Моделирование объемных тел. Моделирование объемных тел слож­ной формы как правило, существенно сложнее, чем плоских, и требует при­менения более сложных элементов. Однако рассмотренный простейший двухузловой элемент без существенных изменений может быть превращен из плоского в объемный. На рис. 2.5 видно, что если плоский элемент со­стоит из двух призм, высота которых равна толщине моделируемого плос­кого тела, то объемный — из двух пирамид. Как призмы, так и пирамиды примыкают с двух сторон к поверхности s, разделяющей соседние ячейки, а ребра призм и пирамид соединяют точки контура поверхности s с узлами 7 и 2. Отрезок 7—2 длиной / перпендикулярен плоскости s, и точка пересечения с этой плоскостью делит его пополам.

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

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

Ячейки, на которые разбита модель, имеют прямоугольную форму, если заданные точки (узлы) 1, 2, 3... расположены ровными рядами, как на рис. 2.3. При произвольном расположении узлов на плоскости или в про­странстве построение ячеек, границы которых равноудалены от ближайших узлов, также возможно, хотя ячейки будут иметь более сложную форму (рис. 2.6).

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

Г раничные условия.

Чтобы рассчитать токи и по­тенциалы по модели, изобра­женной в виде принципиаль­ной схемы на рис. 2.4, необхо­димы граничные условия, т. е. параметры недостающей части электрической цепи, присое­диненной к модели в узлах А и В (см. рис. 2.3 и 2.4).

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

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

Условие второго рода: задан ток или плотность тока на границе. Это соответствует наличию во внешней цепи другого источника питания, под­держивающего постоянный ток, не зависящий от сопротивления модели. В этом случае граница является узлом, потенциал в котором зависит от потен­циала ближайшего внутреннего узла и определяется в ходе решения задачи. Частным случаем такого граничного условия является нулевая плотность тока (адиабатическая граница). Для задач электропроводности это соответ­ствует границе проводника с изолятором. На рис. 2.3 такие условия дейст­вуют на всех внешних границах модели, кроме граничных узлов АиВ.

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

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

• если уравнения моделируемых физических явлений отличаются от линейных;

• если входящие в уравнения коэффициенты не являются константа­ми, а зависят от результатов решения (например, если удельное сопротивле­ние проводника зависит от плотности тока).

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

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

Уравнения теплопроводности и диффузии практически совпадают с уравнениями электропроводности (вместо закона Ома в них используются аналогичные законы Фурье и Фика). Однако скорость протекания этих про­цессов значительно ниже. Вследствие этого не всегда можно ограничиться расчетом установившегося поля температур или распределения примесей, во многих случаях требуется определять распределение параметров во вре­мя переходного процесса. В этом случае решение зависит не только от гра­ничных условий, но и от начального состояния и времени. Соответствующее уравнение отличается от уравнения Лапласа добавлением производной по времени и называется уравнением параболического типа:

82U d2U 82U 8U

ау—— + ау—т + а7—----------- • (2-3)

* дх } 8у2 г 8z 8t К

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

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

Из сопоставления этих двух задач ясно, что для решения нестацио­нарной задачи нужно к первым двум этапам (необходимым и для стацио­нарной задачи) добавить еще два. Таким образом, имеем следующие этапы:

1) уравнение потоков массы или энергии (закон Ома, Фурье, Фика и т. д.);

2) граничные условия 1, 2 и 3-го рода;

3) начальные условия — исходное распределение температур или примесей;

4) условия накопления массы или энергии.

Применительно к тепловой задаче п. 4 соответствует уравнение теп­лоемкости, связывающее теплосодержание с температурой.

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

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

Если шаг решения недостаточно мелкий, то предположение о сохра­нении постоянных тепловых потоков становится некорректным и приводит к неправильному решению (ломаная 16 удаляется от кривой). Максимально допустимый шаг во времени т пропорционален квадрату линейного размера самого мелкого элемента / и зависит от свойств материала (для тепловой задачи он обратно пропорционален коэффициенту температуропроводности а):

Iі

х = —. а

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

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

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

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

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

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

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

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

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