1612725465-eb5831e2a489b993bf8c6c33f84310db (828847), страница 17
Текст из файла (страница 17)
Для этогоаппроксимируем задачу (7.3) разностной схемойµ¶1xj1 +1,0 − xj1 ,0xj1 ,0 − xj1 −1,0w(xj1 +1/2 , 0)− w(xj1 −1/2 , 0)= 0,h1h1h1j1 = 1, . . . , N1 − 1,(7.4)x0,0 = 0, xN1 ,0 = lx .110В результате получилась нелинейная разностная задача для вычисления абсцисс узлов на нижней стороне прямоугольника Ω (ординатывсех узлов равны нулю). Здесь xj1 +1/2 = (xj1 ,0 + xj1 +1,0 )/2. Полученная конечно-разностная задача решается итерационным методом, например, методом последовательных приближений [16]. Аналогично находятся абсциссы xj1 , N2 узлов неравномерной сетки на верхней сторонепрямоугольника Ω, при этом yj1 , N2 ≡ ly .Для поиска координат узлов на левой стороне прямоугольника Ωнеобходимо решить задачуµ¶ddy(0, η)w(0, y)= 0, η ∈ (0, 1),(7.5)dηdηy(0, 0) = 0, y(0, 1) = ly .Конечно-разностный аналог этой задачи записывается какµ¶1y0,j2 +1 − y0,j2y0,j2 − y0,j2 −1w(0, yj2 +1/2 )− w(0, yj2 −1/2 )= 0,h2h2h2j2 = 1, .
. . , N2 − 1,(7.6)y0,0 = 0, y0,N2 = ly ,где yj2 +1/2 = (y0,j2 + y0,j2 +1 )/2. Аналогично определяются ординатыузлов на правой границе (их абсциссы равны xN1 ,j2 ≡ lx ).В результате будут найдены все граничные узлы xj1 ,0 , xj1 ,N2 , x0,j2 ,xN1 ,j2 сетки Ω̄h . Совокупность этих граничных узлов будем обозначатьчерез ∂ Ω̄h . Отметим, что если xj ∈ ∂ Ω̄h , то xj = x(ξj ), где ξ j ∈ ∂ Ξ̄h .7.3. После того как на границе неравномерная сетка построена, вычисляются координаты внутренних узлов адаптивной сетки.
Выведемуравнения метода эквираспределения для определения координат этихузлов. Во-первых, перепишем уравнение Пуассона в новых переменныхξ, η (см. задачу 7.2):µ¶µ¶∂ g22 ∂u g12 ∂u∂g12 ∂u g11 ∂u−+−++ Jf = 0, ξ ∈ Ξ,∂ξJ ∂ξJ ∂η∂ηJ ∂ξJ ∂η(7.7)где J — якобиан преобразования (7.2):J = xξ yη − xη yξ ,g11 = x2ξ + yξ2 ,g22 = x2η + yη2 ,111g12 = xξ xη + yξ yη .Поскольку линейные функции u = x и u = y удовлетворяют уравнению Лапласа, то для произвольного гладкого отображения (7.2) имеемдва тождестваµ¶µ¶∂ g22 ∂x g12 ∂x∂g12 ∂x g11 ∂x−+−+= 0, ξ ∈ Ξ, (7.8)∂ξJ ∂ξJ ∂η∂ηJ ∂ξJ ∂η∂∂ξµg22 ∂y g12 ∂y−J ∂ξJ ∂η¶+∂∂ηµ¶g12 ∂y g11 ∂y−+= 0,J ∂ξJ ∂ηξ ∈ Ξ.(7.9)Далее будем предполагать, что система координат, определяемая искомым отображением (7.2), является ортогональной. Запишем это условие в аналитической форме.
Возьмем произвольные числа ξ0 , η0 из единичного отрезка. При изменении координаты ξ функция x = x(ξ, η0 )будет описывать координатную линию первого семейства с касательным вектором τ 1 = (xξ , yξ ). Координатная линия второго семействаx = x(ξ0 , η), проходящая через точку x(ξ0 , η0 ), имеет касательный вектор τ 2 = (xη , yη ). Поэтому в произвольной точке x(ξ0 , η0 ) условие ортогональности координатных линий τ 1 · τ 2 = 0 можно записать так:g12 (ξ) ≡ 0,ξ ∈ Ξ.(7.10)Кроме того, будем считать, что отображение (7.2) удовлетворяет принципу эквираспределения в дифференциальной форме, аналогичному принципу (2.7.44) в одномерном случае:w(x(ξ))J(ξ) = C = const,ξ ∈ Ξ.(7.11)Использование этих условий в тождествах (7.8), (7.9) приводит к двумерным уравнениям метода эквираспределения:µ¶µ¶∂∂x∂∂xwg22+wg11= 0, ξ ∈ Ξ,∂ξ∂ξ∂η∂η∂∂ξилиµ¶µ¶∂y∂∂ywg22+wg11= 0,∂ξ∂η∂η∂∂ξµ¶µ¶∂x∂∂xk11+k22= 0,∂ξ∂η∂ηгде k11 = wg22 , k22 = wg11 .112ξ ∈ Ξ,ξ ∈ Ξ,(7.12)Итак, в двумерном случае уравнения для определения отображения(7.2) следуют из принципа эквираспределения (7.11) при дополнительном предположении об ортогональности системы координат, задаваемойэтим отображением.Учитывая, что на границе области Ω сетка уже построена, выпишем разностную задачу Дирихле для определения сеточных векторфункций xj — координат внутренних узлов:Λxj = 0, ξ j ∈ Ξh ,xj = x(ξ j ), ξj ∈ ∂ Ξ̄h ,(7.13)где Λ = Λ1 + Λ2 , а разностные операторы Λ1 и Λ2 аппроксимируют совторым порядком соответственно первый и второй дифференциальныеоператоры в левой части уравнения (7.12):µ1 (k11 )j1 +1,j2 + (k11 )j1 ,j2 xj1 +1,j2 − xj1 ,j2·−Λ1 xj =h12h1¶(k11 )j1 ,j2 + (k11 )j1 −1,j2 xj1 ,j2 − xj1 −1,j2−·,2h1µ1 (k22 )j1 ,j2+1 + (k22 )j1 ,j2 xj1 ,j2 +1 − xj1 ,j2Λ2 xj =·−h22h2¶(k22 )j1 ,j2 + (k22 )j1 ,j2−1 xj1 ,j2 − xj1 ,j2 −1−·.2h2Для решения разностной задачи (7.13) можно использовать один израссмотренных ранее итерационных методов, например метод переменных направлений (6.30), (6.31):n+1/2xj− xnjn+1/2= Λ1 xj+ Λ2 xnj ;τ /2(7.14)n+1/2xn+1− xjjτ /2n+1/2= Λ1 xj+ Λ2 xn+1.j(7.15)При этом не надо забывать, что задача (7.13), в отличие от (5.36), является нелинейной, поскольку коэффициенты k11 и k22 операторов Λ1 и Λ2сами зависят от решения xj .
Если эти коэффициенты брать с предыдущей итерации, то реализация дробных шагов (7.14), (7.15) ничем неотличается от линейного случая — используются продольно-поперечныепрогонки.1137.4. Теперь выпишем аппроксимацию задачи (7.1) на криволинейнойсетке Ω̄h . Согласно общей методике, описанной в § 2.7, сначала в исходной задаче надо перейти к новым независимым переменным ξ,η, а затемполученную задачу аппроксимировать на равномерной прямоугольнойсетке Ξ̄h .Задача Дирихле для уравнения Пуассона (7.7) в переменных ξ,η имеет следующий вид:(k11 uξ + k12 uη )ξ + (k21 uξ + k22 uη )η + Jf = 0,u(ξ) = µ( x(ξ) ), ξ ∈ ∂Ξ,гдеk11 =g22,Jk12 = k21 = −g12,Jk22 =ξ ∈ Ξ,(7.16)g11.JВидим, что в новых координатах уравнение Пуассона имеет болеесложный вид, чем в исходных: оно имеет переменные коэффициентыи смешанные производные.
Воспользуемся одной из аппроксимаций этого уравнения, приведенных в работе [12]:Λuj + Jj fj = 0,ξ j ∈ Ξh ,(7.17)где Λ = Λ1 + Λ12 + Λ2 + Λ21 ,1Λ1 uj =h1Λ2 uj =1h2µµk11 (3) + k11 (0) u(3) − u(0) k11 (0) + k11 (1) u(0) − u(1)·−·2h12h1k22 (4) + k22 (0) u(4) − u(0) k22 (0) + k22 (2) u(0) − u(2)·−·2h22h2Λ12 uj =12h1µ¶u(7) − u(6)u(8) − u(5)k12 (3) ·− k12 (1) ·,2h22h2Λ21 uj =12h2µ¶u(7) − u(8)u(6) − u(5)k21 (4) ·− k21 (2) ·2h12h1¶,¶и для нумерации узлов использованы обозначения рис.
2.Для решения разностной задачиΛuj + Jj fj = 0,uj = µ( x(ξ j ) ),114ξ j ∈ Ξh ,ξ j ∈ ∂Ξh(7.18),48CDj27130BA526j1Рис. 2. Шаблон разностного уравненияможно использовать итерационные методы, основанные на схеме приближенной факторизации для уравнений со смешанными производными [12], либо другие, более простые в реализации (но медленнее сходящиеся) итерационные методы, например, метод последовательной верхней релаксации [5]. Чтобы воспользоваться последним методом, перепишем разностные уравнения (7.17) в следующем виде:Ã 8!Xαk uk+ Pj = 0,(7.19)k=0jгде Pj = h1 h2 Jj fj , uk — значение сеточной функции u в узле 9-точечногошаблона, имеющем номер k. Коэффициенты этого уравнения вычисляются по формуламαk =ih2 hk11 (0)+k11 (k) , k = 1, 3,2h1αk =ih1 hk22 (0)+k22 (k) , k = 2, 4,2h2ii1h1hk12 (1) + k12 (2) , α7 = k12 (3) + k12 (4) ,448iiX1h1hαk .α6 = − k12 (2) + k12 (3) , α8 = − k12 (1) + k12 (4) , α0 = −44α5 =k=1Тогда расчетные формулы метода последовательной верхней релаксации примут вид³´.XXnu0 = − Pj +αk un+1+αuα0 ,(7.20)k kkk=1,2,5,6k=3,4,7,8115un+1= τ u0 + (1 − τ )un0 ,0(7.21)где n — номер итерации, τ — итерационный параметр, τ ≥ 1.7.5.
Предположим, что мы хотим сгущать сетку в подобластях больших значений градиента |∇u|. Тогда целесообразно задать управляющую функцию в виде (2.7.53):w(x, y) = 1 + a|∇u|β .(7.22)Однако точное решение u нам не известно, следовательно, невозможновычислить значения управляющей функции и построить адаптивнуюсетку. Поэтому для многомерных задач поступают так же, как в одномерном случае, описанном в § 2.7: необходимо одновременно вестипоиск численного решения uj и построение сетки xj с помощью какойлибо итерационной процедуры.
В качестве начального итерационногоприближения берется равномерная прямоугольная сетка в Ω и на этойсетке решается задача (7.18). Полученное решение uj используется длявычисления управляющей функции в нужных точках. После этого определяется новое положение граничных узлов и путем решения разностной задачи (7.13) строится неравномерная сетка внутри области. Затемна этой сетке вновь решается задача (7.18). Итерационный процесс продолжается до сходимости решения с заданной точностью. Проведенныерасчеты тестовых задач показывают, что зачастую хватает двух – трехитераций, поскольку дальнейшее итерирование уже не приводит к заметному повышению точности численного решения.ЗАДАЧИ7.1. Показать, что если w ≡ 1 на γ, то сетка на границе областибудет равномерной.7.2. Показать, что при преобразовании координат (7.2) уравнениеПуассона (7.1) переходит в уравнение (7.7).7.3.