B.C. Зарубин, д-р техн. наук проф.,

Г.Н. Кувыркин, д-р техн. наук, МГТУ им. Н.Э. Баумана

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

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

Применение современных информационных технологий к проектированию и созданию различных технических объектов нередко приводит к необходимости количественного анализа процессов, протекающих в системах с распределенными параметрами (в континуальных системах) [1]. Информация о распределении в пространстве и изменении во времени температуры, перемещений и деформаций, механических напряжений, скорости и давления жидкости или газа, электрического потенциала, напряженности электрического или магнитного поля и других параметров важна при разработке и оптимизации технологических процессов и рабочих процессов в энергетических установках, при анализе процессов деформирования и динамики конструкций и процессов взаимодействия среды с электромагнитными полями в приборных устройствах. Эта информация может быть получена путем вычислительного эксперимента с использованием математических моделей таких процессов.

Эти модели описываются, как правило, дифференциальными уравнениями с частными производными, решение которых может быть проведено методом конечных разностей (МКР) [2]. Но наряду с МКР для количественного анализа таких моделей перспективен метод граничных элементов (МГЭ) [3] или его сочетание с методом конечных элементов (МКЭ) [4]. Для применения МГЭ математическую модель процесса необходимо предварительно привести к форме, содержащей граничные интегральные уравнения с неизвестными распределениями искомых параметров на границе области, в которой протекает рассматриваемый процесс [1]. Такая модификация математической модели позволяет понизить размерность задачи и тем самым дает возможность сэкономить вычислительные ресурсы. Однако в случае нелинейных процессов неизвестные величины обычно входят не только в интегралы по границе области, но и в интеграл по самой области. В таком случае МГЭ целесообразно сочетать с процедурой последовательных приближений, задаваясь ожидаемым распределением искомых параметров в области и последовательно его уточняя по найденным распределениям этих параметров на ее границе.

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

                                                            (1)

где f(М) - заданная функция положения точки  - дифференциальные операторы Гамильтона и Лапласа. Например, если в однородном теле с постоянным коэффициентом теплопроводности λ произведение λf(М) задает распределение объемной мощности энерговыделения, то уравнение (1) описывает установившееся распределение в этом теле температуры u(М). Если же в однородной среде с постоянной диэлектрической проницаемостью ε произведение εf(М) характеризует объемную плотность электрических зарядов, то уравнению (1) будет удовлетворять распределение в области V потенциала u(М) электрического поля. Этими примерами далеко не исчерпываются приложения (1) к инженерным задачам.

Функция u(М), удовлетворяющая уравнению (1), будет удовлетворять и интегральному соотношению

                               (2)

Функцию w(M₀, M) выберем так, чтобы в случае пространственной (трехмерной) задачи она удовлетворяла уравнению

                                            (3)

где δ(M₀, М) - дельта-функция Дирака, равная нулю в любой точке М, не совпадающей с точкой М₀. При совпадении этих точек дельта-функция неограниченно возрастает, так что для непрерывной функции f(M)

                           (4)

причем Ω(M₀) = 4π, если , если точка М₀ лежит на гладком участке границы F; Ω(M₀) = 0, если  (точка М₀ лежит за пределами области V и ее границы F). Если точка М₀ совпадает с угловой точкой границы F области V, то Ω(M₀) - телесный угол с вершиной в точке М₀, под которым "видна" из точки М₀ внутренность области V.

В случае плоской (двумерной) задачи в уравнении (3) вместо 4π будет множитель 2π, а в уравнении (4) Ω(M₀) = 2π, если , если точка М₀ лежит на гладком участке контура F; Ω(M₀) = 0, если , и, наконец, Ω(M₀) равна внутреннему углу между касательными к контуру в точке М₀, если она является угловой точкой контура.

Пусть r(М₀, М) - расстояние между точками М₀ и М. Тогда функции w(M₀, М) = 1/r(М₀, М) (в случае трехмерной задачи) и w(M₀, М) = -1nr(М₀, М) (в случае двумерной задачи) будут удовлетворять уравнению (3), описывая, например, потенциал электрического поля заряда 1/ε или температурное поле источника тепла мощностью 1/λ, помещенных в точке М₀.

Используя вторую формулу Грина и соотношения (2)-(4), получим граничное интегральное уравнение

                           (5)

где  - единичный вектор внешней нормали к границе F в точке . Правую часть уравнения (5) можно вычислить, если известны распределения u(Р)и q(P) на границей области V. Но корректная постановка задачи может содержать в виде граничных условий в каждой точке  информацию либо о значениях u(Р), либо о значениях q(Р), либо о комбинации их значений, т.е.

                               (6)

Таким образом, исходная информация о граничных значениях u(P) и q(P) недостаточна для непосредственного использования уравнения (5).

Недостающие граничные значения u(Р) и q(P) можно определить приближенно. Для этого на границе области следует выделить N узлов , аппроксимировать распределения u(Р) и q(P) через узловые значения u_n = u(P_n) и q_n = q(P_n) и решить относительно неизвестных N узловых значений систему из N линейных алгебраических уравнений (СЛАУ), которую можно получить из (5) и (6), если в (5) вместо М₀ последовательно полагать Р_n^.. Эта процедура и составляет существо МГЭ (точнее, его прямого варианта), причем (в отличие от МКЭ и МКР) матрица СЛАУ для МГЭ является полностью заполненной.

После вычисления всех узловых значений u_n и q_n на границе F из уравнения (5) нетрудно численным интегрированием найти значение u(М₀) в любой внутренней точке  области V. Численное интегрирование удобно проводить с использованием процедур МКЭ, представив область V системой конечных элементов.

В случае осесимметричной задачи вместо области V достаточно рассматривать половину ее осевого сечения, а вместо границы F - контур этого сечения в координат­ной плоскости 0rz, где r и z - радиальная и осевая координаты цилиндрической системы координат. Тогда задача становится двумерной, а в уравнении (5) w(M₀, М) = 4K(k)/r_*(М₀ М), где K(k) - полный эллиптический интеграл первого рода с модулем k² = 4r(М₀)/r(М)/r²(М₀, М) и r²_*(М₀, М) = (r(М₀).+ r(M))² + (z(M_q) - z(M))². Ясно, что в этом случае число N узловых точек на контуре, а значит, и порядок СЛАУ могут быть существенно уменьшены.

Если рассматриваемый процесс является нестационарным, то (1) переходит в уравнение вида

                                           (7)

где точкой над символом обозначена производная по времени t, а параметр а характеризует скорость выравнивания распределения u(М, t) в области при отсутствии внешних воздействий. Для перехода к математической модели нестационарного процесса в форме, содержащей граничное интегральное уравнение, можно использовать функцию источника [4], однако более гибким является подход, связанный с предварительным переходом в уравнении (7) к конечным разностям по времени. Тогда для момента времени t_k в конце k-го интервала  получим

где u_k(M) = u(M, t_k) и f_k(M) = f(M, t_k) + u(M, t_k-1)/(aΔt_k) а в (5) для момента времени t_k надо w(M₀, M) заменить в случае трехмерной задачи на

и в случае двумерной задачи на

где символом K₀ обозначена модифицированная функция Бесселя второго рода нулевого порядка. Отметим, что на первом интервале времени (k = 1) функция u₀(M) является заданным начальным распределением искомой функции в области V.

Наконец, рассмотрим путь использования МГЭ для моделирования нелинейного нестационарного процесса в неоднородной среде, описываемого уравнением

               (8)

в котором параметры среды λ и c, характеризующие ее свойства проводимости и емкости, зависят не только от положения точки , но и от искомой функции u(M, t).

Аппроксимируя в уравнении (8) производную по времени конечными разностями, в конце k-го интервала времени  получим

где

Тогда для момента времени t_k будет справедливо уравнение (5), но в интеграл по области V вместо f(М) будет входить функция h_k(M, u_k, (M)), зависящая от искомого распределения u_k(М) = u(М, t_k), что, казалось бы, лишает смысла применение МГЭ в случае нелинейного процесса.

Однако сочетание процедуры МГЭ с последовательными приближениями на каждом интервале времени позволяет использовать положительные стороны МГЭ, связанные с понижением размерности задачи. Более того, благодаря полностью заполненной матрице СЛАУ для определения неизвестных граничных значений u_k(P) и , влияние изменения значений искомых величин распространяется на каждой итерации сразу на всю область, что при условии сходимости последовательных приближений существенно сокращает общее число итераций.

Во всех рассмотренных случаях существенным моментом для применения МГЭ была возможность в представляющем математическую модель дифференциальном уравнении выделить в явном виде дифференциальный оператор Лапласа. Это позволяет даже при моделировании нелинейного нестационарного процесса построить достаточно простой алгоритм последовательных приближений, базирующийся на использовании удовлетворяющей уравнению (3) функции w(M₀, M).

Отметим, что в общем случае искомая функция, описывающая процесс в континуальной системе, может быть векторной или тензорной, как, например, при математическом моделировании напряженно-деформированного состояния конструкций [5]. Тогда для применения МГЭ исходную математическую модель необходимо преобразовать к форме, содержащей систему граничных интегральных уравнений, число которых должно совпадать с числом координатных функций, представляющих искомую векторную или тензорную функцию. В этом случае функция w(M₀, M) также будет векторной или тензорной, удовлетворяющей, как правило, уравнению с более сложным, чем в (3) дифференциальным оператором [4], [5].

Список литературы

1. Зарубин B.C. Математическое моделирование процессов в континуальных системах // Информационные технологии. 1995. N 0. С. 11-14.

2. Зарубин B.C., Селиванов В.В. Вариационные и численные методы механики сплошной среды. М.: Изд-во МГТУ, 1993. 360 с.

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

4. Зарубин B.C. Прикладные задачи термопрочности элементов конструкций. М.: Машиностроение, 1985. 292 с.

5. Крауч С, Старфилд А. Методы граничных элементов в механике твердого тела / Пер. с англ. М.: Мир, 1987. 328 с.

ИНФОРМАЦИОННЫЕ ТЕХНОЛОГИИ, №3. 1997

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ