Использование методов подобия и анализа размерностей в компьютерном моделировании

Реферат

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

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

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

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

Математическое моделирование, то есть описание объектов или явлений реального мира с помощью математической символики, является важнейшим видом моделирования. Одной из разновидностей математического моделирования является компьютерное моделирование. Компьютерная модель — это программная реализация математической модели, дополненная различными служебными программами (например, модулями графической визуализации результатов расчетов). Компьютерная модель имеет две составляющие — программную и аппаратную. Программа также является абстрактной знаковой моделью. Такая модель интерпретируется не только людьми, но и компьютером.

Один из путей реализации методологии моделирования является следующий. Первый шаг — это формулировка математической модели ситуации. Такой шаг включает в себя определение относящихся к делу величин и формулировку системы определяющих уравнений, которые описывают проблему в подробностях. На втором шаге мы можем рассматривать и решать задачу, представленную этими модельными уравнениями, как чисто математическую задачу. Когда решение получено, третий шаг состоит в возвращении к исходной постановке и проверки соответствия полученных результатов экспериментальным наблюдениям в исходной практической ситуации. Если имеется соответствие, и решение может помочь в предсказании похожих результатов, мы делаем вывод, что сформулированные модельные уравнения фактически представляют реалистичную модель.

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

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

Методология компьютерного моделирования охватывает различные области — от разработки технических систем до анализа экономических процессов. Вера в то, что прикладная математика и компьютерное моделирование состоят только из разработки методов и алгоритмов для решения задач, возникающих в прикладных областях, является ограниченным и фактически неверным взглядом. Прикладная математика имеет дело со всеми тремя указанными выше шагами, а не только с формальным решением задачи, представленным выше на втором шаге. В действительности, важный аспект прикладной математики и компьютерного моделирования состоит в изучении, исследовании и разработке процедур, являющихся полезными в решении таких задач: это включает аналитические и аппроксимационные методы, численный анализ и методы решения уравнений.

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

Анализ размерностей позволяет нам понять размерностные (то есть связанные с длиной, временем и массой и т.д.) связи величин в уравнениях и результирующие следствия этой размерностной однородности. Масштабирование — это метод, который помогает нам понять величину значений членов уравнений, с помощью сравнения значений величин, появляющихся естественным образом в физических ситуациях.

Метод сравнения размерностей физических величин— это один из основных методов анализа их взаимодействия друг с другом. Уравнения модели должны содержать согласованные между собой переменные величины, то есть быть однородны в смысле размерностей входящих в них величин. Это простое соображение является основой дисциплины, известной как анализ размерностей (Баренблатт).

Основой анализа размерностей является так называемая -теорема. Грубо говоря, эта теорема утверждает, что если имеется некоторый физический закон, задающий связь определенного числа размерных физических величин, то существует эквивалентный закон, который может быть выражен через определенное количество безразмерных величин (традиционно обозначаемых 1, 2, …). Перейдем к более точным формулировкам. (Кутателадзе, Баренблатт)

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

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

Физические величины выражаются числами, которые получаются путем измерения — прямого или косвенного сравнения с соответствующими единицами измерения. Единицы измерения разделяются на основные и производные. Основные единицы измерения задаются в виде тех или иных эталонов, производные единицы измерения получаются из основных в силу определения физической величины.

Фундаментальными физическими мерами являются масса, расстояние, время. При рассмотрении конкретных объектов обычно достаточно этой констатации и предположения о том, что для каждой из этих трех величин можно подобрать некоторые устойчивые значения, которые могут служить в качестве их масштабов, воспроизводимых природными или искусственными эталонами. Такие масштабы массы М, длины L и времени Т называются основными единицами физических величин. (Кутателадзе) Производные единицы конструируются из основных единиц по формуле размерностей

Смысл введения производных единиц заключается в том, что, будучи более тесно связанными с конкретным физическим явлением, они наглядно и компактно характеризуют те или иные достаточно общие физические ситуации. При этом отношение двух численных значений величины не зависит от выбора масштабов. Показатели степеней в формуле для [Zi] могут иметь положительные и отрицательные значения, быть целыми и дробными числами. Однако в силу самого определения размерности эти показатели все одновременно обращаться в нуль не могут.

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

[Z] = L0M0T0 = 1.

(Кутателадзе)

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

кг/M, м/L, с/T,

где M, L, T —положительные числа, показывающие во сколько раз уменьшатся основные единицы массы, длины, времени при переходе от исходной системы СИ к другой системе данного класса. Этот класс обозначается MLT.

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

Размерность физической величины всегда представляет собой степенной одночлен. Этот факт следует из фундаментального общефизического принципа ковариантности (Баренблатт).

Таким образом, следует различать:

а) величины размерные, численное значение которых зависит от принятой системы единиц;

б) величины безразмерные, или отвлеченные, численное значение которых не зависит от принятой системы единиц, [Z] = 1.

Безразмерные величины могут образовываться в виде отношения двух величин одинаковой размерности А

S = Z1/Z2; dim S = А/А = 1

и в виде особой комбинации величин разной размерности в соответствии с формулами

Здесь П — традиционный символ произвольного безразмерного комплекса, а — стандартный символ произведения. Таким образом, комплекс П представляет собой произведение множества величин Zi, имеющих размерности Аi, возведенные в такие степени bi что результирующая размерность оказывается равной 1. Очевидно, что для выполнения этого условия выбор числа размерных величин n не является произвольным, а обусловлен множеством размерностей . Содержательность комплекса П определяется тем, что входящие в него величины отражают конкретные физические взаимодействия. Таким образом, безразмерные комплексы:

а) позволяют уменьшить число переменных, характеризующих изучаемое явление, в том смысле, что заранее устанавливается определенная взаимосвязь между размерными величинами в виде их безразмерных комбинаций, и делают это изучение независимым от системы мер;

б) могут служить мерой, определяющей возникновение критических ситуаций, которые приводят к коренным перестройкам протекания изучаемого явления. [Кутателадзе]

Закономерности, определяемые в физической теории или в эксперименте, всегда можно представить в виде

где величины носят название определяющих параметров, причем аргументы имеют независимые размерности, а размерности аргументов выражаются через размерности определяющих параметров :

Величины имеют независимую размерность, если размерность ни одной из этих величин нельзя представить в виде произведения степеней размерностей остальных величин. Размерность определяемой величины должна выражаться через размерности определяющих параметров :

Если бы это было не так, то размерности величин были бы независимыми и, согласно предыдущему, можно было бы, меняя систему единиц измерения внутри данного класса, произвольно менять величину , оставляя неизменными величины (а следовательно, и все определяющие параметры ). Это означало бы, что величина зависит не только от параметров , т. е. что список определяющих параметров в зависимости заведомо неполон. Таким образом, существуют такие числа , что имеет место формула . Поэтому полагают

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

Следовательно, зависимость представляется на самом деле через функцию аргументов:

или, что то же самое, функция имеет специальный вид

Из этих формул следует, что число безразмерных комплексов, составленных из назначенного набора размерных величин, определяется формулой Бэкингема (Баренблатт)

p=n-k,

где n — полное число переменных; k — число первичных (независимых) размерностей. Этот факт составляет содержание центрального утверждения анализа размерностей — П-теоремы Бэкингема. Безразмерные величины называются параметрами подобия.

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

Пример 1. Физический закон

Связывает расстояние, пройденное материальной точкой при падении в постоянном гравитационном поле, с величиной времени t. В некоторой выбранной системе единиц измерения пусть расстояние x измеряется в см, время t в секундах, а g в см/сек2. Если мы изменяем единицы измерения для основных величин x и t на метры и минуты, то в новой системе единиц

где . Следовательно, g=(длина)(время)-2 и мы получаем

Тогда .

Следовательно, рассматриваемый физический закон не зависит от выбранной системы единиц измерения.

Пример 2. Рассмотрим модель, описывающую падение материальной точки с учетом сопротивления среды. Будем использовать, как обычно, уравнение в форме второго закона механики Ньютона

Напомним, что , где c — безразмерный коэффициент, [S]=L2, []=ML-3. Следовательно, [k]=ML-1.

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

Перепишем исходное уравнение с указанием размерностей слагаемых

Перепишем исходное соотношение в виде

Пусть . Здесь V, — безразмерные переменные, vc, tc — масштабирующие множители соответствующей размерности. Тогда

Преобразуем

Убедимся, что получившиеся в правой части комплексы безразмерны:

Выберем значения vc, tc так, чтобы обеспечить соотношения . Этим соотношениям удовлетворяют следующие значения , .Окончательно, в новых безразмерных переменных V, уравнение принимает вид

Решение дифференциальных уравнений движения в пакете программ MATLAB

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

Найдем решение полученного в безразмерном виде дифференциального уравнения для скорости падающей материальной точки с учетом сопротивления среды на интервале времени [0;1]. После получения решения ответы в размерных единицах можно получить, с помощью преобразований, обратных к проведенным при получении безразмерных переменных.

Для решения задачи составим файл undim.m, который сохраним в рабочей директории пакета. Файл содержит следующее описание функции

function dydt = undim(t,y)

dvdt=1-v^2;

В командном окне зададим следующие команды

>> v0=0;

Это задание начальной скорости.

>> ts=0:.1:1;

>> [T,V]=ode15s(@undim,ts,v0)

Первая из этих команд — задание сетки времени, а вторая — обращение к программе численного решения дифференциального уравнения. В результате получаем

T =

0.1000

0.2000

0.3000

0.4000

0.5000

0.6000

0.7000

0.8000

0.9000

1.0000

V =

0.0995

0.1968

0.2906

0.3792

0.4615

0.5365

0.6040

0.6638

0.7161

0.7615

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

>> plot(T,V)

>> hold on; gtext('график скорости');

Литература

Могилев А.В., Пак Н.И., Хеннер Е.К. Информатика. Учебник для ВУЗов – М.: Издательский центр «Академия», 2007.– 848с.
Лапчик М.П., Семакин И.Г., Хеннер Е.К. Методика преподавания информатики: Учеб. пособие для студ. пед. вузов – М.: Издательский центр «Академия», 2001.-– 624с.
Бенькович Е.С., Колесов Ю.Б., Сениченков Ю.Б. Практическое моделирование динамических систем — СПб.: БХВ-Петербург, 2002. — 464с.
Баренблатт Г.И. Автомодельные явления — анализ размерностей и скейлинг: Учебное пособие — Долгопрудный: Издательский Дом «Интеллект», 2009. — 216 с.
Кутателадзе С. С. Анализ подобия и физические модели.— Новосибирск: Наука, 1986.– 297c.

Движение небесных тел

По закону всемирного тяготения любые два тела, находящиеся на расстоянии r друг от друга и имеющие массы m и M, притягиваются с силой

где G — универсальная гравитационная постоянная. Основываясь на этом законе, опишем движение планет, считая, что m — масса планеты, а M — масса Солнца. Влияние других планет при этом учитывать не будем.

Пусть Солнце находится в начале координатной системы Oxy, а планета в момент времени t находится в точке с координатами x,y. Сила притяжения F, действующая на планету, раскладывается на две составляющие: параллельную оси Ox, равную , и параллельную оси Oy, равную . Используя второй закон Ньютона и его же закон всемирного тяготения, получаем уравнения

Принимая во внимание, что , а , эти уравнения можно переписать в виде

Наконец, учитывая, что , приходим к дифференциальным уравнениям

Запишем эту систему двух дифференциальных уравнений второго порядка в виде системы четырех дифференциальных уравнений первого порядка. Для этого обозначим

Дифференцируя эти соотношения и заменяя соответствующие вторые производные координат по времени, согласно приведенным выше уравнениям Ньютона, получаем

Объединяя полученные результаты, записываем систему в следующем виде

В этой задаче особенно неудобно работать с размерными величинами, измеряемыми миллионами километров, секундами и так далее, поэтому нужно выбрать удобные масштабы измерения величин. В качестве масштаба измерения расстояния удобно брать расстояние от Земли до Солнца =1,4961011м. Чтобы определить удобные масштабы времени и скорости, вычислим параметры движения по круговой орбите с постоянной скоростью. В этом случае сила тяготения выступает в роли центростремительной, которая при постоянной скорости выражается известной из начального курса физики формулой . Таким образом, сила тяготения приравнивается центростремительной силе

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

А при выбранном в качестве r расстоянии от Земли до Солнца и сделанных допущениях, получаем

Таким образом, выбираем в качестве новых переменных следующие величины

Проведем масштабирование (обезразмеривание) уравнений. Для этого подставим в систему уравнений представленные связи новых и старых переменных в каждое из четырех уравнений

, следовательно . Преобразуем получившийся комплекс . Таким образом,

Аналогично,

Далее, , следовательно,

. Преобразуем получившееся уравнение

Следовательно,

Аналогично получается уравнение для Vy. Окончательно получаем систему уравнений

При численном моделировании придется принять некоторое положение условно за начало, а затем изучать движение дальше.

Законы Кеплера

В 1609 году И.Кеплер опубликовал первые два своих закона движения планет. Первый закон описывает форму орбит планет Солнечной системы.

Всякая планета движется по эллиптической орбите с общим фокусом, в котором находится Солнце.

Второй закон характеризует площадь сегмента, заметаемого радиусом–вектором, направленным от Солнца к планете, за некоторые интервалы времени.

За равные промежутки времени радиус–вектор планеты заметает равные площади, так что секториальная скорость постоянна.

Третий закон И.Кеплера, опубликованный в 1619 году, связывает движение различных планет.

Отношение кубов главных осей орбит двух любых планет Солнечной системы равно отношению квадратов периодов их обращения вокруг Солнца.

Напомним, что формула площади криволинейного сектора, граница которого задана в полярных координатах, имеет следующий вид:

Здесь имеется в виду сектор, ограниченный полярной кривой r=f(), и лучами =a, =b.

Получим некоторые равенства, которые понадобятся в дальнейшем. Будем описывать движение планеты с помощью параметризованной кривой

здесь — длина радиуса-вектора, направленного от Солнца к планете, а — угол поворота радиуса-вектора от выбранного фиксированного направления. Будет удобно записывать это в виде

где — двумерный вектор единичной длины соответствующего направления. Отметим, что производная вектора имеет вид и также является вектором единичной длины, который перпендикулярен к . Производная вектора по t — это вектор, имеющий ту же размерность; координаты вектора-производной являются производными координат исходного вектора.

Продифференцируем функцию по правилу дифференцирования сложной функции

и вычислим определитель матрицы, составленный из векторов c(t) и

Получили равенство

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

Дифференцируем это соотношение по t

Объединяя это соотношение с выражением для определителя матрицы, составленного из векторов c(t) и , получим

Поскольку — это скорость заметания радиусом-вектором площади соответствующей области, то второй закон Кеплера эквивалентен утверждению, что , что в свою очередь эквивалентно равенству , но

по свойству определителя. Объединяя полученные результаты, формулируем теорему.

Теорема.

Второй закон Кеплера справедлив тогда и только тогда, когда сила, действующая на планету, является центральной, В этом случае траектория каждой планеты удовлетворяет уравнению

Доказательство.1) Пусть сила, действующая на планету, является центральной. Это может быть выражено равенством (по закону Ньютона сила сонаправлена вектору ускорения, что по свойству определителя дает требуемое). Но последнее равенство эквивалентно равенству , которое, в свою очередь, соответствует равенству

— удвоенной секториальной скорости, то есть справедлив второй закон Кеплера.

2) Пусть выполняется условие постоянства секториальной скорости

Это равенство, однако, можно переписать в виде .

Тогда, дифференцируя левую часть, имеем , а это и есть условие центральности силы. Теорема доказана.

Литература

Могилев А.В., Пак Н.И., Хеннер Е.К. Информатика. Учебник для ВУЗов – М.: Издательский центр «Академия», 2007.– 848с.
Баренблатт Г.И. Автомодельные явления — анализ размерностей и скейлинг: Учебное пособие — Долгопрудный: Издательский Дом «Интеллект», 2009. — 216 с.
Spivak M. Calculus. Publish or Perish, Inc., 1994. – 684p.

Моделирование теплопроводности

Теплота — это форма энергии, которая присутствует в любом веществе. Как любая форма энергии, теплота измеряется в джоулях (1 дж=1 н•м), кроме того, теплота измеряется в калориях ( 1 кал=4,184 дж). Количество теплоты внутри заданного объема определяется с точностью до аддитивной постоянной. Мы будем полагать, что количество теплоты равно нулю, когда температура равна нулю. Предположим, что V — малый объем, внутри которого температура почти одинакова. Экспериментально может быть установлено, что количество теплоты Q пропорционально температуре, которую будем обозначать через u, и массе m=V, где — плотность вещества. Таким образом, количество теплоты задается формулой

Q= u V

Константа c называется теплоемкостью. Она измеряет количество теплоты, требуемое для увеличения температуры на 1 градус в единице объема данного вещества.

Рассмотрим длинный тонкий стержень с теплоизолированной боковой поверхностью. Длина стержня равна L. Зададим координатную ось, которую направим вдоль стержня, координату обозначим через x, при этом 0xL. Так как боковая поверхность стержня теплоизолирована,то отток тепла через нее отсутствует, теплота входит и выходит лишь через концы стержня. Мы предполагаем, что температура u зависти от x в момент t, u(x,t).

Рассмотрим малый участок стержня между x и x+x. Пусть S — площадь поперечного сечения стержня. Объем рассматриваемого участка стержня — Sx, таким образом, количество теплоты равно

Q= cu Sx.

Поэтому количество теплоты в стержне в момент t задается интегралом

Теплоемкость c и плотность вещества иногда меняются от точки к точке. Мы будем, как правило, иметь дело с однородными веществами, в которых обе эти величины являются постоянными.

Уравнение теплопроводности моделирует поток теплоты через вещество. Это уравнение выводится с помощью вычисления изменения количества теплоты Q(t) двумя разными способами. Первый способ — это дифференцирование величины Q(t). Дифференцирование под знаком интеграла дает

Конечно, если теплоемкость и плотность вещества не зависят от времени,

Второй способ вычисления изменения количества теплоты состоит в учете притока и оттока тепла на рассматриваемом участке стержня. Это количество теплоты может измениться за счет притока и оттока тепла через границы участка, а также за источников тепла внутри рассматриваемого участка стержня. Рассмотрим участок стержня между x и x+x. Экспериментальное изучение теплового потока через подобный участок привело к открытию следующих свойств:

тепловой поток направлен от точки с большей температурой к точке с меньшей температурой и пропорционален разности температур ;
тепловой поток пропорционален площади поперечного сечения границ участка;
для теплового потока выполняется закон Фурье: [поток тепла, направленный слева направо]=.

Коэффициент k называется коэффициентом температуропроводности. Поясним смысл выбора знака. Представим, что речь идет о левом торце. Если значение частной производной отрицательно, это означает, что справа от торца стержня температура ниже, чем на торце, и поток будет направлен (от нагретого к менее нагретому) слева направо. Аналогично, если значение частной производной положительно, справа от торца стержня температура выше, чем на торце. Таким образом, поток будет направлен справа налево, а слева направо будет формально направлен отрицательный поток.

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

Здесь f(y,t) плотность источников тепла внутри стержня. Эти источники могут быть образованы, например, за счет проводника тока, расположенного внутри стержня.

Слагаемые внутри квадратных скобок могут быть преобразованы с помощью формулы Ньютона-Лейбница

Если коэффициент k не зависит от x, то

Объединяя полученные формулы, выписываем соотношение, связывающее изменение количества теплоты внутри стержня

Это соотношение выполняется дляx>0 , когда выражение под интегралом равно нулю, то есть

Разделим уравнение на c, обозначим . Окончательно,

Это уравнение теплопроводности, 2 — коэффициент теплопроводности.

Начальные и граничные условия

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

1. (на границе задана температура)

2. (задана температура среды, окружающей торцы)

3. (задан тепловой поток через торцы стержня).

Первый случай как технически, так и идейно, прост и выражается обычно в форме

При фиксированных значениях температуры на концах стержня , . Рассмотрим граничные условия второго типа.

Предположим, что один или оба конца стержня приведены в соприкосновение со средой (жидкостью), температура которой известна. Пусть, для определенности, температура на левом конце равна , а на правом — . Задавая граничные условия этого типа, мы не можем считать значения температуры концов стержня такими же, как у среды, в которую они погружены. Но мы знаем закон Ньютона: если температура одного из концов стержня меньше, чем температура соответствующей жидкости, тепло будет втекать в стержень со скоростью, пропорциональной разности температур. Другими словами

Вытекающий поток тепла при x=0 равен

Вытекающий поток тепла при x=L равен

Здесь h — коэффициент теплообмена. Эти положения вместе с известным законом теплопроводности Фурье можно использовать для получения граничных условий. Закон Фурье был сформулирован выше, для концов стержня все изложенное выглядит так

Вытекающий поток тепла при x=0 равен

Вытекающий поток тепла при x=L равен

Приравнивая выписанные соотношения, получим

0<t<

Обозначая , получаем

Рассмотрим граничные условия третьего типа. Если торцы стержня теплоизоилорованы, то поток через них отсутствует. Это формально записывается так

0<t<

В случае неизолированных торцов стержня задаются потоки, тепла как функции времени

t>0

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

Рассмотрим кратко запись граничных условий в задаче о распределении температуры в тонком теплоизолированном кольце. В этом случае полагаем, что L=2R, где R — радиус кольца. Граничные условия реализуют следующие два равенства: 1)температура при x=0 равна температуре при x=2R; 2)поток, втекающий в «торец» x=0, равен потоку вытекающему из «торца» x=L. Второе из этих соотношений записывается на основе закона Фурье.

Численное решение уравнения теплопроводности

методом конечных разностей

Область решения на плоскости двух переменных для уравнения теплопроводности преобразуется в дискретную сетку из узлов Здесь — целые числа. Узлы сетки ; , где — шаги сетки по координатам и соответственно, — целые числа.

Неизвестная функция , участвующая в краевой задаче, заменяется сеточной функцией на узлах сетки. Частные производные по координатам заменяются соответствующими конечными разностями. Пусть — шаг сетки вдоль рассматриваемой координаты, — значение функции в рассматриваемой точке , а и последующее и предыдущее значения функции по данной координате. Тогда первая производная по этой координате может быть заменена следующей конечной разностью

а вторая производная выглядит следующим образом:

Для первой производной по времени конечная разность имеет вид

А для второй производной по координате конечная разность выглядит следующим образом:

Чтобы написать разностную схему, приближенно описывающую решаемую задачу, нужно выполнить следующие два шага.

1. Необходимо заменить область непрерывного значения аргумента областью его дискретного изменения.

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

После осуществления такой процедуры мы приходим к алгебраической системе уравнений. Таким образом, задача о численном решении исходной начально-краевой задачи для дифференциального уравнения теплопроводности сводится к вопросу о нахождении решения полученной алгебраической системы. Остановимся на этих вопросах подробнее.

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

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

Пример 1. Равномерная сетка на отрезке. Разобьем единичный отрезок [0,1] на равных частей. Расстояние между соседними узлами назовем шагом сетки. Точки деления — узлы сетки.

Пример 2. Равномерная сетка на плоскости. Рассмотрим множество функций двух аргументов . В качестве области определения выберем прямоугольник

Разобьем отрезки [0,1] оси и [0,T] оси соответственно на и частей, пусть . Через точки деления проведем прямые, параллельные соответствующим осям. В результате пересечения этих прямых получим узлы , которые и образуют сетку. Эта сетка имеет шаги h и соответственно по направлениям x и t.

Область изменения переменных задачи следующая {0tT, 0xL}. Заменяем непрерывные переменные дискретными:

— начальные условия,

— граничные условия на левом конце, — граничные условия на правом конце.

Заменяем частные производные их конечно-разностными аппроксимациями

Подставим эти выражения в уравнение и разрешим получившееся соотношение относительно ui+1,j, то есть значения функции в верхнем временном слое

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

Граничные условия имеют следующие сеточные аналоги. Условиям теплоизолированных концов

очевидно, должны соответствовать условия вида

ui+1,1=ui+1,2, ui+1,N=ui+1,N-1

Условия, соответствующие заданным потокам,

заменяются, соответственно, для левого конца

тогда .

Аналогично, для правого конца

Откуда .

При рассмотрении задачи о распределении температуры в тонком теплоизолированном кольце граничные условия в конечно-разностном виде приводят к системе двух линейных уравнений с двумя неизвестными. При этом первое уравнение выглядит так

ui+1,N=ui+1,1

Значение, которому равны эти величины, легко вычисляется.

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

Следовательно, .

Окончательно, .

Для правого торца выражение определяется аналогичным образом.

Численное моделирование теплопроводности в пакете MATLAB

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

Рассмотрим решение этого уравнения на интервале времени 0t4. Переменная x изменяется в пределах -5x5. Граничные условия представлены в виде на левом конце, на правом конце стержня. Начальные условия для распределения температуры были заданы в виде функции

Чтобы решить представленную задачу, нужно создать m-файл следующего вида

function u=heat(alpha,x,t,init,bdry)

%решение одномерного уравнения теплопроводности

%векторы x и t с начальным условием u(x,t(1))=init

%и граничными условиями

%u(x(1),t)=bdry(1),u(x(end),t)=bdry(2)

N=length(x);

M=length(t);

dx=mean(diff(x));

dt=mean(diff(t));

s=alpha*dt/dx^2;

u=zeros(M,N);

u(1,:)=init;

for i=1:M-1

u(i+1,2:N-1)=s*(u(i,3:N)+u(i,1:N-2))+...

(1-2*s)*u(i,2:N-1);

u(i+1,1)=bdry(1);

u(i+1,N)=bdry(2);

end

который следует сохранить в рабочей директории под именем heat.m.

После этого в командном окне пакета набираются следующие команды:

>>tvals=linspace(0,4,81);

>>xvals=linspace(-5,5,21);

>>init=20+5*sign(xvals);

>>uvals=heat(2,xvals,tvals,init,[15 25]);

>>surf(xvals,tvals,uvals)

>>xlabel x; ylabel t; zlabel u

Первые две команды формируют дискретную сетку значений на заданных интервалах изменения переменных задачи, третья команда формирует массив начальных условий. Далее следует обращение к функции, сохраненной в m-файле, а последние две команды обеспечивают визуализацию полученного решения (см. рисунок).

Литература

1. Несис Е.И. Методы математической физики. Учебн. пособие для студентов физ.-мат. фак. пед. ин-тов. М., «Просвещение», 2007.

2.Дьяченко В.Ф. Основные понятия вычислительной математики. М., «Наука», 2012.

3.Данилина Н.И., Дубровская Н.С., Кваша О.П., Смирнов Г.Л., Феклисов Г.И. Численные методы. Учебник для техникумов. М., «Высшая школа», 2006. 368с.

4. Hunt Brian R. Matlab R2007 с нуля®! Книга + Видеокурс.: / — М.: Лучшие книги, 2008. — 352 с.:

Использование методов подобия и анализа размерностей в компьютерном моделировании