Кирьянов Д. - MathCad 11 (1077323), страница 53
Текст из файла (страница 53)
13.6. Для дискретизации второй производной по пространственной координате необходимо использовать три последовательных узла, вто время как для разностной записи первой производной по времени достаточно двух узлов. Записывая на основании данного шаблона дискретноепредставление для (i,k)-ro узла, получим разностное уравнение:= Dk+1-11-Ci-1Рис. 13.6. Шаблон аппроксимации явной схемыдля уравнения теплопроводностиПриведем в разностной схеме (8) подобные слагаемые, перенеся в правуючасть значения сеточной функции с индексом к (как часто говорят, с предыдущего слоя по времени), а в левую — с индексом k+i (т. е. со следующего временного слоя).
Кроме этого, введем коэффициент с , который будетхарактеризовать отношение шагов разностной схемы по времени и пространству_ 2 • At • р 1 Л . Несколько забегая вперед, заметим, что значение параметра с, называемого коэффициентом Куранта, имеет большое значение для анализа устойчивости разностной схемы. С учетом этихзамечаний, разностная схема (8) запишется в виде:С •ui,k+i-— ui-i,k/\+ I1 - C i . k K kС-+ — ui+i,k+i,k»(9)Множители для каждого из значений сеточной функции в узлах шаблона,соответствующие разностному уравнению (9), приведены рядом с каждойточкой шаблона на рис. 13.6.
Фактически, геометрия шаблона и эти множители задают построенную нами разностную схему.326Часть III. Численные методыНесложно убедиться в том, что для получения замкнутой системы разностных алгебраических уравнений систему (9) необходимо дополнить дискретным представлением начального и граничных условий (6) и (7). Тогда числонеизвестных будет в точности равно числу уравнений, и процесс формирования разностной схемы будет окончательно завершен.ПримечаниеВажно подчеркнуть, что возможная нелинейность полученной системы алгебраических уравнений определяется зависимостями от температуры функцийD(u) и ф(и), т. е.
как коэффициент диффузии, так и источник тепла могут бытьфункциями сеточной функции u i k .Если присмотреться к разностным уравнениям (9) повнимательнее, то можно сразу предложить несложный алгоритм реализации этой разностной схемы. Действительно, каждое неизвестное значение сеточной функции со следующего временного слоя, т. е.
левая часть соотношения (9) явновыражается через три ее значения с предыдущего слоя (правая часть), которые уже известны. Таким образом, в случае уравнения теплопроводностинам очень повезло — для расчета 1-го слоя по времени следует попроступодставить в (9) начальное условие (известные значения и с нулевого слоя вузлах сетки), для расчета 2-го слоя достаточно использовать вычисленныйтаким образом набор и с 1-го слоя и т. д. Из-за того, что разностная схемасводится к такой явной подстановке, ее и называют явной, а благодаря пересчету значений с текущего слоя через ранее вычисленные слои - схемойбегущего счета.Линейное уравнениеСделанные замечания относительно реализации явной схемы для уравнениядиффузии тепла сразу определяют алгоритм ее программирования вMathcad.
Для решения задачи нужно аккуратно ввести в листинг соответствующие формулы при помощи элементов программирования.Решение системы разностных уравнений (9) для модели без источников тепла, т. е. ф(х,т,и=о и постоянного коэффициента диффузии D=const приведена в листинге 13.1. В его первых трех строках заданы шаги по временной и пространственной переменным т и Д, а также коэффициентдиффузии D, равный единице. В следующих двух строках заданы начальные(нагретый центр области) и граничные (постоянная температура на краях)условия, соответственно.
Затем приводится возможное программное решение разностной схемы, причем функция пользователя v(t) задает векторраспределения искомой температуры в каждый момент времени (инымисловами, на каждом слое), задаваемый целым числом t.Глава 13. Уравнения в частных производных327\ Листинг 13.1. Явная схема для линейного уравнения теплопроводностиТ := 0 . 0 0 0 5М := 20МD { u ) := 1ф ( х , и ) := 0Init(х):=ФBorder (t)m : - 0 ..(х-0.45)-Ф(х-0.55):= 0Мu m : = I n i t (m -Д)F(v):=v l э <- Border (т • т) + 0 . • {vo + vi)VIM <- B o r d e r (t • т) + 0 . • (v M + v M - l )f o r m e 1 ..M- 1,l | ( A,)D(v m _i)-T+D{VIW.I)-Tvm-vm+ivlV ( t ) :=иift = 0F ( V ( t - 1) )otherwiseРИС. 13.7. Решение линейного уравнения теплопроводности(листинг 13.1)328Часть HI.
Численные методыНачальное распределение температуры вдоль расчетной области и решениедля двух моментов времени показано на рис. 13.7 сплошной, пунктирной иштриховой линиями, соответственно. Физически такое поведение вполнеестественно — с течением времени тепло из более нагретой области перетекает в менее нагретую, а зона изначально высокой температуры остывает иразмывается.Нелинейное уравнениеНамного более интересные решения можно получить для нелинейногоуравнения теплопроводности, например с нелинейным источником теплаФ(и)=103- (u-V). Заметим, что в листинге 13.1 мы предусмотрительно определили коэффициент диффузии и источник тепла в виде пользовательскихфункций, зависящих от аргумента и, т.
е. от температуры (если бы собирались моделировать явную зависимость их от координат, то следует ввести впользовательскую функцию в качестве аргумента переменную х, как этосделано для источника тепла ф). Поэтому нет ничего проще замены определения этих функций с констант D ( U ) = I И ф(х,и)=о на новые функции, которые станут описывать другие модели диффузии тепла. Начнем с того, чтопоменяем четвертую строку листинга 13.1 на ф(х,и) =ю3- (u-u 3 ), не изменяяпока постоянного значения коэффициента диффузии.(~Примечание^С физической точки зрения зависимость коэффициента диффузии и функцииисточника тепла от температуры означает, что эти параметры будут менятьсяот точки к точке среды, определяясь локальными значениями текущей температуры в этих точках. Ввод ненулевого источника тепла означает, что среда получает определенное количество тепла, тем большее, чем больше локальнаятемпература.
Можно догадаться, что введение такой зависимости может моделировать, в частности, горение среды.Если осуществить расчеты с упомянутым источником (имеющим кубическую нелинейность), то получится очень интересное решение уравнениятеплопроводности, имеющее профиль тепловых фронтов. С течением времени граница раздела высокой и низкой температуры распространяется вобе стороны от зоны первичного нафева, оставаясь весьма четко выделенной (рис. 13.8).Еше более неожиданные решения возможны при нелинейности также и коэффициента диффузии. Например, если взять квадратичный коэффициентдиффузии D(X,U}=U (что с учетом его умножения на неизвестные функциисоздаст кубическую нелинейность уравнения), а также ф{х,и) = ю 3 и 3 5 , тоВы сможете наблюдать совсем иной режим горения среды.
В отличие отрассмотренного эффекта распространения тепловых фронтов, горение оказывается локализованным в области первичного нагрева среды, причем,температура в центре нафева со временем возрастает до бесконечной вели2Глава 13. Уравнения в частных производных329чины (рис. 13.9). Такое решение описывает так называемый режим горения«с обострением»., Ч > <jV£,G-G-©0-5VI20UОI I- „0.5Рис. 13.8. Решение уравнения теплопроводностис нелинейным источником (тепловой фронт)Читателю предлагается поэкспериментировать с этим и другими нелинейными вариантами уравнения теплопроводности.
Существенно, что такиеинтересные результаты удается получить лишь численно, а в Mathcad только с применением элементов программирования.111Xi1•>('»1i•ii<*A.\г..т.кГ.Е15..ЭГ..0-4I11*i11*1КiV'(»a \* 1- . « . l A ^ У-.,ЛftЛ ЛЛЛт> • и у * в у "ч*> ^ J *№ ^£^ ^ У " * U~fО.бРис. 13.9. Решение уравнения теплопроводности с нелинейным источникоми коэффициентом диффузии (режим локализации горения)330Часть III.
Численные методыУстойчивостьКак мы убедились, явная разностная схема Эйлера дает вполне разумныерезультаты и вполне может использоваться для практического моделирования задач, связанных с решением уравнений в частных производных. Однако теперь пришло время сказать об очень важной характеристике разностных схем, которая называется их устойчивостью.
Не вдаваясь в детали,заметим, что производить расчеты можно только при помощи устойчивыхразностных схем, а чтобы пояснить это понятие, обратимся вновь к листингу 13.1, реализующему явную схему для линейного уравнения диффузии.Слегка изменим соотношение шагов по времени и пространственной координате, произведя расчеты сначала с т=о.ооо5 (эти результаты показаны нарис. 13.7 выше) и также с т=о.оою и T=o.ooi5. Результат применения разностной схемы Эйлера для т=о.оою примерно тот же, что и для меньшегозначения т., приведенного на рис.
13.7. А вот следующее (казалось бы, незначительное) увеличение шага по времени приводит к катастрофе(рис. 13.10). Вместо ожидаемого решения получается совершенно неожиданные профили температуры, которые быстро осциллируют вдоль пространственной координаты, причем амплитуда и число пиков этих осцилляции быстро увеличиваются от шага к шагу. Совершенно ясно, чтополученное решение не имеет ничего общего с физикой моделируемого явления, а является следствием внутренних свойств самой разностной схемы,которые до этого были для нас скрыты.ф-ч0.5Рис.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
