XIII Власова Б.А., Зарубин B.C., Кувыркин Г.Н. Приближенные методы математической физики (1081417), страница 81
Текст из файла (страница 81)
Обращением матрицы С систему (11.33) можно привести к нормальной системе ОДУ и для ее решения использовать один из вариантов метода Рунге — Кутты. Однако обращение матрицы С далеко не всегда оправдано с точки зрения экономии вычислительных ресурсов ЭВМ при полном решении задачи, особенно в тех случаях, когда существенна зависимость теплоемкостн тела от температуры и эту зависимость необходимо учесть путем решения (11.33) итерациями, последовательно уточняя элементы матрицы С . Рассмотрим некоторые способы приближенного решения системы (11.33), имея в виду для нелинейных задач возможность учета зависимости от температуры не только теплоемкости, но и коэффициента теплопроводности тела. Используем конечно-разностную аппроксимацию производной в (11,33) в пределах интервала времени Ь1л =1л — 1я 1, приняв д  — В д1 Л1ь где индексы. Й вЂ” 1 и Й соответствуют моментам времени в начале и конце к-го интервала.
Применим двухслойную разностную схему с весами и вместо (11.33) получим (11.34) где ц б 10, 1]. Индексы /с — 1 и к в обозначениях матриц С и Л указывают на то, что элементы этих матриц в общем случае 1П2. Задачи телдовроводности а твердом теде 589 е„"= е„, + (с„,) '(д», — л е„,)л1 . В правой части этого равенства все параметры известны или могут быть вычислены по найденным на предыдущем итервале времени (при Й = 1 — по начальным) узловым значениям температуры, составляющим матрицу-столбец Е» . Следует иметь в виду, что использование явной схемы накладывает ограничение на выбор Ь1» иэ условия устойчивости (см.
8.3). При и= 1 (11.34) принимает вид е* е„= ( —," +л») (с;,— '+о;). (11.35) Выбирая ~)И», можно руководствоваться лишь соображениями точности расчета, поскольку неявнал схема устойчива при любых значениях Ы». В случае существенной зависимости теплоемкости и коэффициента теплопроводности тела от температуры вычисления по (11.35) на каждом интервале времени приходится проводить несколько раз, последовательно уточняя значения элементов матриц С» и Л».
При слабой зависимости этих элементов от температуры обычно достаточно ограни- не остаются постоянными. В частности, они могут изменяться при изменении температуры в силу возможной зависимости от нее теплоемкости и коэффициента теплопроводности тела. В предельных случаях О = О и О = 1 (11.34) соответствует явной и неявной разностным схемам аппроксимации (11.33). Но в отличие от метода конечных разностей (МКР) выбор ц= О не приводит к распаду (11.34) на независимые уравнения относительно неизвестных узловых значений температуры.
Сохраняется лишь преимущество явной схемы, связанное с возможностью решения нелинейных задач без итераций на каждом интервале времени, поскольку после обращения матрицы С», можно разрешить (11.34) явно относительно искомой матрицы- столбца 590 Пс ПРИКЛАДНЫЕ ЗАДА ЧИ читься лишь первым приближением, приняв Е.„= ( — '-'+Л„,) '(С„., ',-'+д;). (И.зб) Отметим, что в случае нелинейных граничных условий и зависимости объемной мощности 1г энерговыделения от темпера- И) туры элементы матрицы-столбца Я также заранее неизвестны при расчете на Й-и интервале н возникает необходимость их последовательного уточнения в (11.35) и (11.36).
Если в (11.24) приближенно принять вт(1,М) т,(м) — т„,(м) д1 Ь|ь рассматривал правую часть (11.24) в момент времени 1ь, то в соответствии с методом прямых получим последовательность краевых задач для дифференциально-разностного уравнения су(л(м)сут,(м)) - — т,(м)+9,(м) =0, м ~1, й ~ 1ч, с(м) Ыл с граничными условиями Ть(Р) = ~1(1ю Р), Р Е ом и 'л(Р) гутя(Р) п(Р) + Ц(Р)ть(Р) = Д,(Р), Р Е Яз — — Я'1 Яы где оь(м) =ф1(1ь,м)+ ( ) " '( ), ~ь(Р) =Явь,Р) и Бь(Р) = Ь4 = ~з(1ь, Р). Каждой иэ этих задач соответствует функционал (ХЦ у(т,)=~ ~'-(чт,) + — т„-д„т„~~ и +~ ~~ — т„-у,т,)ия, ,/ ~2 минимизируемый на распределениях температуры Ть(М), удовлетворяющих граничному условию Ть(Р) = ~1(1л, Р), Р Е Яд, //.2. Эадачи теалоароводности в твердом теле 591 непрерывных на замыкании 1/ = 1" О о' и имеющих кусочно непрерывные производные в области т', ограниченной поверхностью 5. Используя процедуру МКЭ, этот функционал можно привести к виду (11.18) н затем из условий его минимума получить СЛАУ вида (11.20), решение которой будет эквивалентно (11.35).
Но в случае нелинейной задачи, минимизируя такой функционал, например, методом локальных вариаций, можно непосредственно найти приближенное распределение Т»(М) (М Е т') температуры в момент времени 1». Выбор // = 1/'2 в (11.34) приводит к двухслойной сил/д/е/прачкой разносганой схеме и более высокому порядку погрешности по Ь/» по сравнению с остальными значениями //. Того же эффекта можно добиться, если для решения задачи нестационарной теплопроводности испольэовать КЭ в пространственновременнбй области", однако и в том, и в другом случае необходимо уточнять элементы матриц С, Л" и Я при решении нелинейных задач.
Избежать последовательного уточнения элементов этих матриц можно с помощью неявной трехслойной разнесенной схемы сс" -6" тт +6 +6 » /с — 2 +Л /с» 1 /с 2 д (11 37) 2Ь/ ' "' 3 устойчивой при любых значениях Ы. Здесь зависящие от температуры элементы матриц С», Л» и Я», вычисляют по известным на Й-м интервале времени элементам матрицы-столбца с/» /, соответствующим моменту времени /» / в середине удвоенного интервала 2Ы = /» — /» 2. Поэтому (11.37) можно явно разрешить относительно искомой матрицы-столбца ет». Отметим, что при й = 1 из-за отсутствия информации об элементах матрицы-столбца /В» для нахождения /с//» придется использовать одну из двухслойных разностных схем.
'См., например: Зврдонн В.С., Селиванов В.В. 592 11. ПРИКЛАДНЫЕ ЗАДА ЧИ 11.3. Двумерное течение вязкой жидкости др ~~ ~д иг Ь; — — +ц'~ — '=О, дх; ~-~ дхг 1=1 1 г 1=1,2; ~~~ — =О, (11.38) ди; дхг где б, и и( — проекции на оси Ох; векторов плотности и скорости объемных сил соответственно; р — давление; коэффициент сдвиговой вязкости жидкости. Разобъем область г' на К лагранэсевых конечных элементов (КЭ), образующих сетку КЭ с общим числом Жн узлов. На этой сетке введем систему (ун)н функций вида (11.17) е и, р (М) =~) ~ Г)~'~~ р(')(М), (У = 11, 1ч'е, (11.39) ею1 п=1 где чэ„(М), и = 1, Ԅ— функции формы КЭ с номером е = (е) = 1, Е и числом узлов Ж,; Й„н — элементы матрицы Й, размера (в) 11)1 х Жн, устанавливающей соответствие между номерами узлов в отдельном КЭ и в сетке КЭ (см.
10,3). Приняв (11.39) в качестве проекционных функций, в соответствии с условиями (6.76) ортогональности невязки Же-мерному подпространству Ик гильбертова простронства Я непрерывных на 7 функций, имеющих в г' кусочно непрерывные Метод конечных элементов (МКЭ) можно применить к решению задач, математическая формулировка которых включает несколько искомых функций и в нее входит более чем одно дифференциальное уравнение. В качестве примера, связанного с решением системы уравнений, рассмотрим достаточно медленное установившееся двумерное течение однородной несжимаемой вязкой жидкости параллельно координатной плоскости х10хг (см.
3.2) в области г', ограниченной кусочно гладким контуром Г. Такое движение описывает система трех уравнен и й вида (3. 22) 1КЗ, двумерное течение ввзной в(ндкоетн 593 производные, умножим (11.38) на функцию (рь(М), Ь = 11, оп, н проинтегрируем по области Е: (рь(М) ~6((М) — — +т~"7~~()((М)) ИР=О, 1=1,2, (11.40) др(М) х; ( (рь(М) ) ' ЙР=О, А=1,1Уп. (11.41) два(М) р зо( Здесь ((7з — двумерный оператор Лапласа в плоскости х~Охз. Функции (рь(М) имеют в области Р кусочно непрерывные производные по координатам х;, 1= 1, 2.
Поэтому, используя первую у)орлеулу Грина, запишем е!(н)ч ч(м)о =~е!(Р)(чч(Р)) (Р)я— ( ее Г (т7(ру,(М))Чи(М) ЙГ, 1= 1, 2, где (7 н п(Р) — двумерный оператор Гамильтона и единичныв вектор внешней нормали в точках Р Е Г контура Г соответственно. Тогда вместо (11.40) получим (е! (и ) (); (м) — — ) — ц (те, (м) ) о.; (и) ) не + др(М) дх; +п) е!(Р)(о;(Р)) (Р)ИГ=О, !=!,2, !.=!,ш!. (!!.42) Г Распределения давления и проекций скорости в области Е аппраксимируем при помощи их узловых значении рв( и ив); соотношениями р(М) = ",) рог(руч(М), гн(М) = ~ и((г(~р((г(М), 1 = 1, 2, Ы. ПРИКЛАДНЫЕ ЗАДАЧИ 594 которые подставим в 111.41) уь(М)~~ оу;~> с)Р'=О, 1=1,2, 1,=1,%е, х; 1м1 А =! и (11.42) Мх — П ~Ярс(М))им ~~) ~у~дм(М)Ы+ %=1 Мх +ц/рь(Р)~з~;(Р)) (Р)м=О, '=1,2, Е=!,м.
г ч=1 Отсюда следует система линейных алгебраических уравнений 1СЛАУ) в матричной записи (11.43) В11 =Я, ~зА =Рм, 1=1,2, Р1= 1,Же. 1~з1м-ц+, = мяъ Матрица-столбец Я имеет тот же размер и злементы где à — матрица-столбец размера ЗМЕ х 1, злементами которой являются 595 1ЬЗ. двумерное течение ввзкой жидкости и Язс = О, г = 1, 2, Ь = ГУ~. Элементами симметрической матрицы Р порядка ЗАс~ являются 1'~ дру.(М) др,ч(М) Рз1с ц+ьз1м,1+,— — 0б;, ) ~~ с1г, р /ж1 дд„(М) Рз1~.- 1+,зх = Рзь,з(м- 1+ = Ю(М) д с1~' Х1 и Рзп з,ч = О, 1, 1 = 1, 2, Ь, Ас = 1, Ж~, где бо —— 1 и ри 1 = 1' и 4 = 0 при е~1. Отметим, что второй интеграл в выражении для Яз1у.-цаз отличен от нуля, если узел с номером Ь принадлежит участку границы Г, на котором заданы значения (~ун;)в.