Методы компьютерного моделирования сварочных процессов

Глава 12

Методы компьютерного моделирования сварочных процессов

 

Уравнение сохранения заряда. Первое правило Кирхгофа

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

.                                                 (12.4)

Например, для узла 1

;

;

;

,                             (12.5)

где .

Мы получили уравнение (12.5), в котором несколько неизвестных – потенциал узла 1 и смежных с ним узлов 2, 3, 4 и 5. Это - уравнение постоянства заряда в узле 1. Такое уравнение можно записать для каждого узла, поэтому число уравнений равно общему количеству неизвестных. Вместе они образуют систему линейных уравнений:

,                        (12.6)

которая в сокращенной (тензорной) форме записи выглядит как

,                                  (12.7)

где - вектор неизвестных потенциалов в узлах,  - матрица системы уравнений, - порядок матрицы (число неизвестных и уравнений в системе). При тензорной записи пропускают знак суммы. На то, что в выражении (12.7) записан не одночлен, а сумма одночленов, указывает «слепой» - повторяющийся два раза индекс k, по которому и производят суммирование.

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

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

.

Можно сказать, что каждое уравнение в системе (12.6) – это проинтегрированное по некоторой области (по одной ячейке) дифференциальное уравнение (12.2).

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

Свойства системы уравнений

Матрица традиционно называется матрицей жесткости (потому что МКЭ создавали для задач механики), но в данном случае это матрица проводимости . Эта матрица симметричная:

.

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

Действительно, если потенциалы смежных ячеек равны

,                                               (12.8)

то все токи будут равны нулю. Подставим (12.8) в уравнение (12.5)

,

откуда

,

.

Таким образом, матрица является положительно определенной. Это обеспечивает устойчивость решения системы методом Гаусса.

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

Очень важным моментом является однородность системы уравнений (в правой части - нули). Это значит, что она всегда имеет решение . Кроме того, мы установили, что  также будет решением при любом . Такие решения называются тривиальными.

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

Если в другом узле задать, что , получим разные потенциалы во всех узлах.

Условием 2-го рода является подача тока от внешнего источника в какой-то из узлов. Этот ток добавляется к правой части соответствующего уравнения.

 

Устойчивость решения

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

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

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

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

Рис. 12.4. Схема роста погрешности решения при вырождении матрицы

 

Итерационные методы решения

Возможны два пути применения итераций при решении системы уравнений.

I. Итерационное решение системы уравнений (применяется вместо метода Гаусса). Порядок итераций такой (см. схему на рис. 12.3):

1) Задаемся начальными значениями потенциала в точках 1, 2, 3 …, например  (кроме А, В,…, в которых потенциал  известен).

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

- у точки 1 в уравнение постоянства заряда входят значения в точках со второй по пятую; поскольку все они имеют нулевой потенциал, значение в точке 1 после первой итерации не изменится: ;

- у точки 2 среди соседей есть точки A и B, потенциалы которых  и  известны, а , по этим данным из уравнения  постоянства заряда можно найти ;

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

3) Теперь будем считать  новыми значениями в т. 1, 2, 3 и повторим итерационное уточнение значений потенциала:

- в точке 1 потенциал изменится , поскольку на предыдущей итерации изменился потенциал в смежной с ней точке 2;

- в точке 2 потенциал не изменится , поскольку на предыдущей итерации не изменялся потенциал в смежных с ней точках A, B, 1 и 3;

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

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

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

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

II. Итерационное уточнение решения, полученного методом Гаусса. Чтобы организовать итерационный процесс, нужно представить искомое значение на итерации в виде суммы значения на предыдущей итерации и приращения за итерацию:

.                                       (12.12)

Подставив это выражение в систему уравнений (12.7), получим:

;

.                                    (12.13)

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

Ход итерационного процесса такой:

1) Задаемся начальными значениями потенциала .

2) Находим из решения системы (12.13) , затем по формуле (12.12).

3) Если после подстановки этих значений в исходные уравнения (12.7) получаем тождества, то найденные значения являются точными, и итерации прекращаются. Если же невязка  ненулевая, то выполняется следующая итерация. Невязка от предыдущей итерации оказывается в правой части уравнений (12.13), и мы ищем такую добавку к значениям неизвестных , чтобы она скомпенсировала невязку от предыдущей итерации. Суммируя предыдущее решение с этой добавкой, находим новое приближение по формуле (12.12).

4) Условием завершения итераций является уменьшение до заданного порога либо невязки , либо приращения за шаг .

Уравнение Фурье

Уравнение Фурье

,                                         (12.16)

является аналогом закона Ома (12.1): вместо плотности электрического тока  в нем присутствует удельная мощность теплового потока , вместо электропроводности  - коэффициент теплопроводности , а вместо напряженности (градиента) электрического поля  – градиент температур на границе  ( – вектор нормали к границе).

Теплота распространяется от горячих частей тела к холодным (в сторону убывания температуры), поэтому в уравнении присутствует знак минус. Поток тепла пропорционален градиенту температур, т. е. разности температур на единицу длины в направлении потока теплоты, и зависит от свойств материала (теплопроводности). Размерность теплопроводности – Вт/(м·К).

Если имеется ячейка (элементарный объем) и найдены потоки теплоты через все ее границы, то общее количество теплоты, поступившего в объем за время

,                                   (12.17)

где   – удельная мощность потока теплоты через границу .

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

.                                               (12.18)

Изменение температуры при изменении теплосодержания зависит от теплоемкости вещества. Удельная теплоемкость – это количество теплоты, необходимое для повышения температуры единицы массы вещества на , размерностью Дж/(кг·К).

;                                                (12.19)

;

.                                      (12.20)

Уравнение теплопроводности (уравнение Пуассона) имеет вид:

,                                   (12.21)

где  – накопление теплоты,  – изменение теплосодержания.

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

- выравнивание температуры между точками тела;

- остывание (выравнивание температуры между телом и окружающей средой).

Температура передается от точки к соседним точкам, поэтому при моделировании применима явная схема: следующее состояние точки можно определить по текущему состоянию данной и соседних с ней точек.

 

Метод элементарных балансов

В литературе по численным методам расчета тепловых процессов существует понятие «модель из теплоемких масс и теплопроводящих стержней». Речь идет именно о тех принципах, которые реализованы в программном комплексе «Сварка». Ту электрическую схему, которой мы заменяли деталь из электропроводного материала, можно назвать «моделью из электропроводящих стержней».

Каждый элемент (граница ячеек = электропроводящий стержень) соединяет две ячейки (два узла). Нашу тепловую модель можно представить такими эквивалентными схемами (рис. 12.11).

Рис. 12.11. Схемы моделей для задач энергомассопереноса: из конденсаторов и сопротивлений (а); из масс и стержней (б):   R - сопротивление, C - электроемкость конденсатора, - теплопроводность, - теплоемкость

 

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

 

Неявная схема решения

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

Согласно (2.24), в каждой внутренней ячейке

,                 (12.26)

где N – количество границ ячейки. В этом выражении, в отличие от (12.24), потоки тепла записаны не для начала, а для конца текущего шага решения. Если ячейка имеет наружные границы с условиями 2 и 3 рода, то добавляются потоки тепла через эти границы:

. (12.27)

Получаем систему уравнений относительно неизвестных температур конца шага

.                                  (12.28)

Диагональные коэффициенты матрицы проводимости

;         (12.29)

внедиагональные коэффициенты матрицы

;                           (12.30)

Правая часть – вектор поступившей теплоты

.         (12.31)

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

.         (12.32)

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

Теория деформаций

Существует три вида движения твердого тела (рис. 12.14): поступательное, вращательное и деформационное (изменение расстояний между точками тела). Изменение линейных размеров называется линейной деформацией, а изменение угловых размеров - угловой.

Рис. 12.14. Схема перемещений и деформаций твердого тела

 

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

;

.

Компоненты линейной деформации выражают через изменения длины отрезков в теле (разности перемещений их концов):

;

 и т. д.

Общая формула:

.                                          (12.29)

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

.

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

[рад];

.

Введя обозначение , получим единую формулу для малых линейных и угловых деформаций, которая включает в себя формулу (12.29), как частный случай:

.                        (12.30)

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

.                                       (12.31)

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

Рис. 12.15. Искажение деформации при повороте тела

 

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

Формула, исключающая ложные деформации от поворота тела или его фрагмента в пространстве:

.                                (12.32)

В этой формуле применена сокращенная тензорная запись скалярного произведения векторов, выражение  записывают как , пропуская знак суммы. О необходимости вычисления и суммирования произведений при всех значениях «слепого» индекса i сигнализирует его двукратное повторение в формуле. Традиционная запись последнего члена уравнения (12.32)

.

Для линейной деформации из (12.32) получаем

.

Формула (12.32) дает компоненты деформации в системе координат, связанной с волокнами материала. В результате перемещений точек эта система поворачивается в пространстве относительно исходной системы координат. Компоненты поворота, найденные по формулам

                                        (12.33)

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

Деформацию, измеренную по отношению к исходным размерам тела , называют деформацией Коши,

,                                             (12.34)

а по отношению к размерам тела после его деформирования  - деформацией Альманси:

,

.                                  (12.35)

Истинные деформации, или деформации Генки, определяют по отношению к текущим размерам тела, меняющимся в процессе деформирования:

.

Для этого необходимо разбить процесс на бесконечно мелкие шаги и найти деформацию интегрированием от  до :

.              (12.36)

Различие между тремя видами деформации становится существенным при .

Деформации Коши и Альманси несимметричны при растяжении и сжатии. Если отрезок длиной  растянуть в два раза, то , , , а . Если тот же отрезок сжать в два раза, то , , , а . Деформация Коши  означает, что размер стал нулевым , а дальше размер меняет знак, «выворачивается наизнанку». Аналогично поведение деформации Альманси при . Деформацию Генки называют истинной, поскольку она имеет одинаковые числовые значения при растяжении и сжатии в два раза и сохраняет физический смысл при любых значениях деформации:  имеет место, когда тело вытягивается в длинную нить, а , когда тело расплющивается в тонкую пластину. Возможность пересчета любого из трех видов деформации в остальные показывает, что все они равноценны.

Совокупность компонент деформации называют тензором:

.                                      (12.37)

Тензор – это вектор второго порядка, характеризуемый величиной и парой направлений (каждая компонента деформации связана с двумя осями, вдоль одной оси направлен отрезок, вдоль другой – перемещение конца отрезка). Закон парности касательных деформаций, являющийся следствием формул (12.31) и (12.32):

                                                 (12.38)

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

Первый инвариант тензора – относительное изменение объема :

,

где  - символ Кронекера, .

Шаровой тензор соответствует состоянию всестороннего равноосного растяжения или сжатия, при котором изменяется объем тела без изменения его формы:

.                                  (12.39)

Найдя первый инвариант любого тензора, можно выделить из него шаровую составляющую (шаровой тензор):

.                      (12.40)

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

;                                        (12.41)

;           (12.42)

.                       (12.43)

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

 при ,

,                               (12.44)

то такие оси и компоненты деформации называют главными:

.                               (12.45)  

У шарового тензора все оси главные.

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

Рис. 12.16. Связь линейных и угловых деформаций

 

Всего тензор имеет шесть независимых параметров, в том числе три инварианта. Полный комплект инвариантов включает три параметра. Примеры таких комплектов: ; . Любые комбинации из инвариантов также являются инвариантами.

Главные деформации  – это корни кубического уравнения

                        (12.46)

или

.                               (12.47)

Способ нахождения главных компонент:

1) имеются формулы расчета инвариантов  по компонентам произвольного тензора, например:

;

2) найдя инварианты, нужно составить и решить кубическое уравнение (12.47).

Один из важнейших инвариантов девиатора деформации – интенсивность деформации:

.            (12.48)

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

Первый инвариант деформации показывает изменение объема эллипсоида (независимо от формы и соотношения осей), второй - суммарное изменение формы (сумму квадратов отклонений точек эллипсоида от эквивалентной по объему сферы). Третий инвариант, например параметр Лодэ

                                             (12.49)

показывает особенности формы эллипсоида: если  (осевое растяжение), то ; если  (кручение), то ; если  (двухосное растяжение), то .

Деформированное состояние можно изобразить в виде кругов Мора (рис. 12.17). Значение первого инварианта (объемной деформации) отображает расположение кругов на оси ; значение второго (интенсивности деформации) – радиус большого круга, значение третьего (параметра Лодэ ) - соотношение размеров малых кругов.

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

 

Теория напряжений

Механическим напряжением σ называется отношение силы P, действующей между двумя частями твердого тела, к площади A разделяющей их поверхности (рис. 12.18):

.                                                  (12.50)

Рис. 12.18. Схема определения нормального и касательного напряжений

 

По размерности напряжение аналогично давлению и измеряется в мегапаскалях: .

Если поверхность плоская, а сила перпендикулярна плоскости (действует по нормали к ней), напряжение называется нормальным:

,

а если сила действует вдоль плоскости (по касательной), напряжение называется касательным:

.

Нормальная сила аналогична силе давления, а касательная – силе трения. В общем случае сила P направлена под углом к поверхности и может быть разложена на составляющие  и . Напряженное состояние в элементарном объеме характеризуют компоненты напряжения на трех секущих поверхностях, перпендикулярных осям координат. Силу, передаваемую через поверхность, можно разложить по осям координат и найти для каждого сечения три компоненты напряжения: нормальную и две касательных. Для обозначения компоненты используют два индекса, один из которых обозначает ось, нормальную к поверхности, а другой – ось, по которой действует составляющая силы. В тензорной записи оси обозначают номерами 1, 2 и 3, нормальные компоненты , а касательные компоненты . Общая формула