Учебное пособие 2049

(распределение тепла, электростатических и магнитостатических полей, безвихревое течение идеальной жидкости и др.). Дискретизация одного из уравнений этого класса – уравнения Пуассона – подробно рассмотрена в разделах II – III.

Типичным примером параболического уравнения является уравнение теплопроводности

где k(x, y, z) – коэффициент теплопроводности, Q(x, y, z, t) – количество тепла, генерируемого в ед. объема, – плотность вещества, (x, y, z, t) –распределение температуры. Это урав-

нение дополняется начальными и краевыми условиями: = 0

при t=t0 ; = на 1 и –k n=q на 2.

Рассмотрим конечно-элементную формулировку нестационарной задачи теплопроводности для стержня 0 x L (одномерный случай). Для простоты положим =с=1. Поскольку такая задача не имеет естественного вариационного принципа, воспользуемся проекционным методом, в частности методом Галеркина.

Потребуем, чтобы проекции невязок на интервале [0, L] и на границе x=L на базисные функции {Ni} и {Wi} были равны нулю.

31

Построим аппроксимацию для : = i(t)Ni(x). Ограничив выбор Wj=–Nj при x=L и Nj=0 при x=0, получим систему

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

В качестве весовых функций возьмем -функции Дирака

где tn tn+1 – tn, – некоторое число. Для неизвестной функции на интервале [tn, tn+1] воспользуемся аппроксимацией

Здесь n (tn), n+1 (tn+1).

Потребовав, чтобы

и положив = 1/2, что соответствует наиболее часто используемой схеме Кранка-Николсона, приходим к матричному уравнению относительно n+1

32

уравнений (1.39) для определения вектора n+1. Решение таких систем осуществляется методом Ньютона-Рафсона, который сводится к итерационному процессу, на каждом шаге которого решается система линейных алгебраических уравнений.

В заключение приведем пример гиперболического уравнения. Рассмотрим случай гибкой струны, на которую действует постоянное натяжение . Концы струны закреплены, и она совершает малые колебания около положения устойчивого равновесия – интервала 0 x 1 оси x.Если (x, t) – перемещение точки струны перпендикулярно x, то для кинетической и потенциальной энергий можно соответственно записать

где – плотность струны и c2 = / . Согласно принципу Гамильтона, интеграл

t1

F( ) (T V)dt

t0

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

22

c2t2 x2

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

33

Упражнения

1. Провести конечно-элементную дискретизацию уравнения Пуассона

( ) = – .

Выделив два соседних элемента, показать, что при этом

условие 1 1 2 2 на их общей границе выполняется ав-

n n

томатически.

2.Как в методе конечных элементов учесть условие на

резкой границе раздела сред, заданное в виде Dn2 – Dn1 = ( – поверхностная плотность заряда) для электростатической задачи?

3.Получить дискретные уравнения задачи теплопроводности для тонкого однородного стержня

(x, 0) = 10, (0, t) = 10 + 0.2t, (a, t) = 10– 0.3t, где (x, t) –

температура стержня в точке x в момент времени t. k принять равным 1.

4. Концы струны закреплены в точках x=0 и x=1; первоначально струна находится в покое ( / t=0 при t=0) и имеет

форму = (x). Найти численное решение, описывающее последующие перемещения струны. Построить дискретные уравнения для МКЭ, пользуясь вариационным и проекционным представлением задачи.

34

2. КОНЕЧНО-ЭЛЕМЕНТНАЯ ДИСКРЕТИЗАЦИЯ УРАВНЕНИЯ ЛАПЛАСА С ДОПОЛНИТЕЛЬНЫМИ ОГРАНИЧЕНИЯМИ

В стандартной краевой задаче для уравнения Лапласа

дополнительно ставятся граничные условия 1-го рода

и/или 2–3 рода

на границе 2 (2.3)

n

(случай 0 соответствует условиям 2-го рода). В эквивалентной вариационной формулировке решение данной задачи соответствует минимуму функционала

причем условие (3.25)(2.2) учитывается как главное, а (2.3) – как естественное. Очевидно, что если граничный интеграл в (2.4) отсутствует, то на границе 2 будет выполняться однородное условие Неймана

0 .n

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

В ряде случаев на потенциал задачи (2.1)-(2.3) могут накладываться дополнительные ограничения [2, 5, 32]. Рассмотрим подробнее построение дискретной конечно-элементной модели

35

с учетом таких ограничений. При этом для простоты исключим из рассмотрения граничные условия 1-3 рода; их учет не влечет никаких особенностей и соответствует стандартной процедуре (см. раздел 1.4).

2.1. Учет скачка потенциала на разрезе

Пусть на некоторой поверхности разреза S ставится условие

что соответствует заданию скачка потенциала на S. Здесь + и– –значения неизвестной функции в одной и той же точке по разные стороны разреза.

Минимизируя функционал задачи

( i – коэффициенты, Ni – базисные функции – полиномы), придем к системе уравнений

Sij i 0, i

где Sij Ni Njd . Матрица {Sij} является симметричной и

положительно определенной.

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

36

Следуя стандартному методу конечных элементов [28], проведем объединение этих систем.

Если S не является линией разреза, то в узлах 2 и 3 должно выполняться условие непрерывности, т.е.

(21) (22) 0,

(31) (32) 0.

В этом случае получается система (с учетом обозначений

i i(l))

Если на S происходит скачок, то каждому узлу, лежащему на S, необходимо поставить в соответствие два значения потенциала i+ и i– . Одно из этих значений будет относиться к элементам, лежащим по одну сторону от S, другое – к элементам, лежащим по другую сторону относительно S. В нашем случае условие скачка можно записать в виде

(1)2 (2)2 I,

(1)3 (2)3 I.

37

При объединении элементов для каждого узла на S следует оставить только одну независимую переменную, например, 2 (21), 3 (31) , а (22) и (32) заменить соответственно

выражениями 2 – I, 3 – I. В итоге получим систему

S11(1) 1 S12(1) 2 S13(1) 3 0,

S21(1) 1 S22(1) S22(2) 2 S23(1) S23(2) 3 S24(2) 4 IS22(2) IS23(2), S31(1) 1 S32(1) S32(2) 2 S33(1) S33(2) 3 S34(2) 4 IS32(2) IS33(2),

S42(2) 2 S43(2) 3 S44(2) 4 IS42(2) IS43(2).

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

4 0.

имеет отличный от нуля определитель. В случае, если I = 0 (скачка нет), ее решение – тождественный нуль, если I 0 (скачок есть), получим искомое ненулевое решение задачи.

38

39

в качестве узлового значения узла 2 следует брать 2, а внутри

элементов 2 и 4 – величину ( 2 – I).

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

Осесимметричная формулировка отличается от плоской только способом вычисления матричных элементов Sij(e) :

Sij(e) r Ni(e) N(je)drdz.

e

2.2. Учет потока градиента потенциала через заданную поверхность разреза

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

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

Получим вариационную формулировку такой задачи. Пусть функционал имеет вид

где I ( + – –)|S – скачок потенциала на поверхности разреза

S; в различных точках этой поверхности значения потенциала могут заметно отличаться, но разность + и –, вычисленных в одной точке, остается постоянной. Поэтому слагаемое Ф0I в

40