Методы вычислений
 Скачать 2.39 Mb.
|
|
w(x, Итак, пусть задана управляющая функция w(x, y) ≥ 1 и требуется построить адаптивную сетку на границе γ области Ω. Рассмотрим, для определенности, нижнюю сторону прямоугольника Ω. В соответствии с методом эквираспределения нам нужно решить задачу (2.7.39): d dξ µ w(x, 0) dx(ξ, 0) dξ ¶ = 0, ξ ∈ (0, 1), x(0, 0) = 0, x(1, 0) = Получить решение x = x(ξ, 0) нелинейной задачи (7.3) в виде аналитической формулы удается очень редко. Поэтому для поиска отображения) применим конечно-разностный метод. Для этого аппроксимируем задачу (7.3) разностной схемой, 0) x j 1 +1,0 − x j 1 ,0 h 1 − w(x j 1 −1/2 , 0) x j 1 ,0 − x j 1 −1,0 h 1 ¶ = 0, j 1 = 1, . . . , N 1 − 1, (7.4) x 0,0 = 0, x N 1 ,0 = l x . 110 В результате получилась нелинейная разностная задача для вычисления абсцисс узлов на нижней стороне прямоугольника Ω (ординаты всех узлов равны нулю. Здесь x j 1 +1/2 = (x j 1 ,0 + x j 1 +1,0 )/2. Полученная конечно-разностная задача решается итерационным методом, например, методом последовательных приближений [16]. Аналогично находятся абсциссы x j 1 , узлов неравномерной сетки на верхней стороне прямоугольника Ω, при этом y j 1 , N 2 ≡ Для поиска координат узлов на левой стороне прямоугольника необходимо решить задачу, y) dy(0, η) dη ¶ = 0, η ∈ (0, 1), y(0, 0) = 0, y(0, 1) = l y . (7.5) Конечно-разностный аналог этой задачи записывается как, y j 2 +1/2 ) y 0,j 2 +1 − y 0,j 2 h 2 − w(0, y j 2 −1/2 ) y 0,j 2 − y 0,j 2 −1 h 2 ¶ = 0, j 2 = 1, . . . , N 2 − 1, (7.6) y 0,0 = 0, y 0,N 2 = где y j 2 +1/2 = (y 0,j 2 + y 0,j 2 +1 )/2. Аналогично определяются ординаты узлов на правой границе (их абсциссы равны x N 1 ,j 2 ≡ В результате будут найдены все граничные узлы x j 1 ,0 , x j 1 ,N 2 , сетки ¯ Ω h . Совокупность этих граничных узлов будем обозначать через ∂ ¯ Ω h . Отметим, что если x j ∈ ∂ ¯ Ω h , то x j = x(ξ j ), где ξ j ∈ ∂ ¯ Ξ h 7.3. После того как на границе неравномерная сетка построена, вычисляются координаты внутренних узлов адаптивной сетки. Выведем уравнения метода эквираспределения для определения координат этих узлов. Во-первых, перепишем уравнение Пуассона в новых переменных, η см. задачу 7.2): ∂ ∂ξ µ g 22 J ∂u ∂ξ − g 12 J ∂u ∂η ¶ + ∂ ∂η µ − g 12 J ∂u ∂ξ + g 11 J ∂u ∂η ¶ + Jf = 0, ξ ∈ где J — якобиан преобразования (7.2): J = x ξ y η − x η y ξ , g 11 = x 2 ξ + y 2 ξ , g 22 = x 2 η + y 2 η , g 12 = x ξ x η + y ξ y η . 111 Поскольку линейные функции u = x и u = y удовлетворяют уравнению Лапласа, то для произвольного гладкого отображения (7.2) имеем два тождества 0, ξ ∈ Ξ, (7.8) ∂ ∂ξ µ g 22 J ∂y ∂ξ − g 12 J ∂y ∂η ¶ + ∂ ∂η µ − g 12 J ∂y ∂ξ + g 11 J ∂y ∂η ¶ = 0, ξ ∈ Далее будем предполагать, что система координат, определяемая искомым отображением (7.2), является ортогональной. Запишем это условие в аналитической форме. Возьмем произвольные числа ξ 0 , из единичного отрезка. При изменении координаты ξ функция x = x(ξ, будет описывать координатную линию первого семейства с касательным вектором τ 1 = (x ξ , y ξ ). Координатная линия второго семейства = x(ξ 0 , η), проходящая через точку x(ξ 0 , η 0 ), имеет касательный вектор. Поэтому в произвольной точке x(ξ 0 , η 0 ) условие ортогональности координатных линий τ 1 · τ 2 = 0 можно записать так) ≡ 0, ξ ∈ Кроме того, будем считать, что отображение (7.2) удовлетворяет принципу эквираспределения в дифференциальной форме, аналогичному принципу) в одномерном случае) = C = const, ξ ∈ Использование этих условий в тождествах (7.8), (7.9) приводит к двумерным уравнениям метода эквираспределения: ∂ ∂ξ µ wg 22 ∂x ∂ξ ¶ + ∂ ∂η µ wg 11 ∂x ∂η ¶ = 0, ξ ∈ Ξ, ∂ ∂ξ µ wg 22 ∂y ∂ξ ¶ + ∂ ∂η µ wg 11 ∂y ∂η ¶ = 0, ξ ∈ или 0, ξ ∈ где k 11 = wg 22 , k 22 = wg 11 112 Итак, в двумерном случае уравнения для определения отображения) следуют из принципа эквираспределения (7.11) при дополнительном предположении об ортогональности системы координат, задаваемой этим отображением. Учитывая, что на границе области Ω сетка уже построена, выпишем разностную задачу Дирихле для определения сеточных вектор- функций x j — координат внутренних узлов 0, ξ j ∈ Ξ h , x j = x(ξ j ), ξ j ∈ ∂ где Λ = Λ 1 + Λ 2 , а разностные операторы и аппроксимируют со вторым порядком соответственно первый и второй дифференциальные операторы в левой части уравнения (7.12): Λ 1 x j = 1 h 1 µ (k 11 ) j 1 +1,j2 + (k 11 ) j 1 ,j2 2 · x j 1 +1,j 2 − x j 1 ,j 2 h 1 − − (k 11 ) j 1 ,j2 + (k 11 ) j 1 −1,j2 2 · x j 1 ,j 2 − x j 1 −1,j 2 h 1 ¶ , Λ 2 x j = 1 h 2 µ (k 22 ) j 1 ,j2+1 + (k 22 ) j 1 ,j2 2 · x j 1 ,j 2 +1 − x j 1 ,j 2 h 2 − − (k 22 ) j 1 ,j2 + (k 22 ) j 1 ,j2−1 2 · x j 1 ,j 2 − Для решения разностной задачи (7.13) можно использовать один из рассмотренных ранее итерационных методов, например метод переменных направлений (6.30), (6.31): x n+1/2 j − x n j τ /2 = Λ 1 x n+1/2 j + Λ 2 x n j ; (7.14) x n+1 j − x n+1/2 j τ /2 = Λ 1 x n+1/2 j + При этом не надо забывать, что задача (7.13), в отличие от (5.36), является нелинейной, поскольку коэффициенты и операторов и сами зависят отрешения. Если эти коэффициенты брать с предыдущей итерации, то реализация дробных шагов (7.14), (7.15) ничем не отличается от линейного случая — используются продольно-поперечные прогонки 7.4. Теперь выпишем аппроксимацию задачи (7.1) на криволинейной сетке ¯ Ω h . Согласно общей методике, описанной в § 2.7, сначала в исходной задаче надо перейти к новым независимым переменным ξ,η, а затем полученную задачу аппроксимировать на равномерной прямоугольной сетке Задача Дирихле для уравнения Пуассона (7.7) в переменных ξ,η имеет следующий вид+ k 12 u η ) ξ + (k 21 u ξ + k 22 u η ) η + Jf = 0, ξ ∈ Ξ, u(ξ) = µ( x(ξ) ), ξ ∈ где k 21 = Видим, что в новых координатах уравнение Пуассона имеет более сложный вид, чем в исходных оно имеет переменные коэффициенты и смешанные производные. Воспользуемся одной из аппроксимаций этого уравнения, приведенных в работе [12]: Λu j + J j f j = 0, ξ j ∈ где Λ = Λ 1 + Λ 12 + Λ 2 + Λ 21 , Λ 1 u j = 1 h 1 µ k 11 (3) + k 11 (0) 2 · u(3) − u(0) h 1 − k 11 (0) + k 11 (1) 2 · u(0) − u(1) h 1 ¶ , Λ 2 u j = 1 h 2 µ k 22 (4) + k 22 (0) 2 · u(4) − u(0) h 2 − k 22 (0) + k 22 (2) 2 · u(0) − u(2) h 2 ¶ , Λ 12 u j = 1 2h 1 µ k 12 (3) · u(7) − u(6) 2h 2 − k 12 (1) · u(8) − u(5) 2h 2 ¶ , Λ 21 u j = 1 2h 2 µ k 21 (4) · u(7) − u(8) 2h 1 − k 21 (2) · u(6) − и для нумерации узлов использованы обозначения рис. Для решения разностной задачи+ J j f j = 0, ξ j ∈ Ξ h , u j = µ( x(ξ j ) ), ξ j ∈ ∂Ξ h (7.18) 114 Рис. 2. Шаблон разностного уравнения можно использовать итерационные методы, основанные на схеме приближенной факторизации для уравнений со смешанными производными, либо другие, более простые в реализации (но медленнее сходящиеся) итерационные методы, например, метод последовательной верхней релаксации [5]. Чтобы воспользоваться последним методом, перепишем разностные уравнения (7.17) в следующем виде+ P j = где P j = h 1 h 2 J j f j , u k — значение сеточной функции u в узле точечного шаблона, имеющем номер k. Коэффициенты этого уравнения вычисляются по формулам 2h 1 h k 11 (0)+k 11 (k) i , k = 1, 3, α k = h 1 2h 2 h k 22 (0)+k 22 (k) i , k = 2, 4, α 5 = 1 4 h k 12 (1) + k 12 (2) i , α 7 = 1 4 h k 12 (3) + k 12 (4) i , α 6 = − 1 4 h k 12 (2) + k 12 (3) i , α 8 = − 1 4 h k 12 (1) + k 12 (4) i , α 0 = Тогда расчетные формулы метода последовательной верхней релаксации примут вид − ³ P j + X k=1,2,5,6 α k u n+1 k + X k=3,4,7,8 α k u n k ´. α 0 , (7.20) 115 u n+1 0 = τ u 0 + (1 − τ где n — номер итерации, τ — итерационный параметр, τ ≥ 1. 7.5. Предположим, что мы хотим сгущать сетку в подобластях больших значений градиента |∇u|. Тогда целесообразно задать управляющую функцию в виде (2.7.53): w(x, y) = 1 + Однако точное решение u нам неизвестно, следовательно, невозможно вычислить значения управляющей функции и построить адаптивную сетку. Поэтому для многомерных задач поступают также, как вод- номерном случае, описанном в § 2.7: необходимо одновременно вести поиск численного решения и построение сетки с помощью какой- либо итерационной процедуры. В качестве начального итерационного приближения берется равномерная прямоугольная сетка вина этой сетке решается задача (7.18). Полученное решение используется для вычисления управляющей функции в нужных точках. После этого определяется новое положение граничных узлов и путем решения разностной задачи (7.13) строится неравномерная сетка внутри области. Затем на этой сетке вновь решается задача (7.18). Итерационный процесс продолжается до сходимости решения с заданной точностью. Проведенные расчеты тестовых задач показывают, что зачастую хватает двух – трех итераций, поскольку дальнейшее итерирование уже не приводит к заметному повышению точности численного решения. З АДА Ч И. Показать, что если w ≡ 1 на γ, то сетка на границе области будет равномерной. Показать, что при преобразовании координат (7.2) уравнение Пуассона (7.1) переходит в уравнение (7.7). 7.3. Показать, что при w ≡ 1 линейное отображение x = ξl x , y = является решением уравнения (7.12). 7.4. Доказать, что при w ≡ 1 координаты узлов равномерной на сетки являются решением задачи (7.13). 7.5. Показать, что разностная схема (7.18) аппроксимирует задачу) со вторым порядком пои. Метод конечных элементов. Метод конечных элементов является эффективным способом численного решения различных задач для дифференциальных уравнений с частными производными со сложной формой области решения. Для упрощения изложения рассмотрим суть этого метода на том же простейшем примере, который использовался в предыдущем параграфе задаче Дирихле для уравнения Пуассона с двумя пространственными переменными в прямоугольнике Ω (7.1). Как ив одномерном случае, вначале сведем эту задачу к задаче с однородными краевыми условиями. Для этого введем функцию v = u − u 0 , взяв в качестве функцию, которая удовлетворяет граничным условиям задачи (7.1). Тогда функция v является решением задачи + ˜ f (x) = 0, x ∈ Ω, v(x) = 0, x ∈ γ = где ∆ — оператор Лапласа (x, y) = f (x, y) + y l y µ xx (x, l y ) + ¡ 1 − y l y ¢ µ xx (x, 0)+ + x l x µ yy (l x , y) + ¡ 1 − x l x ¢ µ yy (0, Таким образом, вместо задачи (7.1) можно решать задачу + f (x) = 0, x ∈ Ω, u(x) = 0, x ∈ Введем линейный оператор A = −∆ с областью определения ∈ C 2 ( ¯ Ω), u ¯ ¯ γ = Тогда задача (8.1) может быть записана в виде операторного уравнения = Используя скалярное произведение в L 2 (Ω) и взяв произвольные функции u, v ∈ D A ⊂ L 2 (Ω), согласно формуле Грина, получаем, v) = − Z Ω v∆u dxdy = Z Ω (u x v x + u y v y ) dxdy = (u Av), 117 поэтому оператор A является симметрическим [16] на D A . Для произвольной функции u ∈ справедливы равенства, u) = Z Ω (u 2 x + u 2 y ) dxdy = Z Ω |∇u| 2 dxdy = поэтому из оценки (см. задачу 8.1) ||u|| L 2 (Ω) ≤ min(l x , l |