Изучение распространения тепла в пластине

 

Кафедра прикладной математики

 

                                                                         

 

 

Учебная дисциплина: Вычислительная математика

 

 

Пояснительная записка  к курсовой работе.

Тема: “Изучение распространения  тепла в пластине”.

Вариант №2.

 

 

 

                      Студент:

                      Преподаватель:

                                                                                                                                                      

                                                                                                                                                                      

 

 

 

                           

Оглавление

Раздел:      стр:

1. Исходные физические и математические модели. Постановка задачи.   3
1. Нелинейная модель с распределёнными параметрами.    3
2. Нелинейная модель с сосредоточенными параметрами.    4
3. Линейная модель с сосредоточенными параметрами.    4
4. Исходные данные.         5
2. Построение оценки зависимости коэффициента конвективного теплообмена от

температуры.           6

3. Определение момента установления температуры окружающей среды.   7

   3.1 Уточнение корня методом  половинного деления.     8

   3.2 Уточнение корня комбинированным  методом.     9

   3.3 Уточнение корня методом итераций.       11

4. Вычисление интеграла I.         12

   4.1 Вычисление интеграла  J по формуле прямоугольников.    12

   4.2 Вычисление интеграла  J по формуле трапеций.     13

   4.3 Вычисление интеграла  J по формуле парабол (Симпсона).    14

5. Приближённое решение задачи Коши.       15

   5.1 Решение задачи Коши  методом Эйлера.      16

   5.2 Решение задачи Коши  методом Рунге-Кута.      17

6. Вывод.            18

 

 

 

1. Исходные физические и математические модели. Постановка задачи.

На практике часто встречается задача исследования температурного поля внутри тел различной формы. Пусть однородная пластина толщины D, имеющая в начальный момент времени температуру U₀, помещена в среду с температурой Q, изменяющейся во времени по заданному закону Q = Q(t). Считая толщину пластины малой по сравнению с остальными размерами пластины, можно рассматривать температуру ее внутренних точек как функцию

U=U(t), где x

[-D/2,D/2] – абсцисса данной точки пластины,

t≥0 – время.

Очевидно, что в силу симметрии задачи U(x,t) – чётная функция аргумента x, поэтому её достаточно рассматривать при x [0,D/2].

1. Нелинейная модель с распределёнными параметрами

 

Функция U(x,t) будет представлять собой решение уравнения теплопроводности:

 

 

 

удовлетворяющими начальному условию                      

 

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

        , где

– коэффициент температуропроводности материала пластины

- коэффициент теплопроводности материала пластины.

- коэффициент конвективного  теплообмена между пластиной  и окружающей среды

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

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

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

Будем считать, что коэффициент конвективного теплообмена зависит от температуры пластины по закону:

, где - известные функции, а

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

 

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

   , где - некоторая функция, а – известные коэффициенты; известные функции

функция , определённая при t≥0, и величина Q₀ — заданы, а момент t₀ достижения установленной температуры определяется как корень уравнения

j(t) = Q₀ ( 2-й этап курсовой работы)

 

 

2. Нелинейная модель с сосредоточенными параметрами.

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

Задача Коши:

 - нелинейное диф. уравнение 1 порядка и начальные условия.

Полученную нелинейную задачу Коши можно решить методами Рунге-Кутта и Эйлера. (Решение этой задачи - 4-й этап курсовой работы).

 

3. Линейная модель с сосредоточенными параметрами.

Будем считать, что

  где - постоянная времени системы.

Решение линейной задачи Коши:

 

Это решение рассмотрим при  :

Обозначим     для

 

 

Приближённое вычисление интеграла I 3-й этап курсовой работы. Сравнивая этот результат с решением задачи Коши, можно оценить влияние, оказываемое на процесс установления температуры переменным характером коэффициента α. Если окажется, что это влияние незначительно и на практике им можно пренебречь, то, очевидно, и при решении исходной задачи будет оправданным считать .

 

4. Исходные данные.

 

A [м2/ч]	λ [ккал/м*ч*°С]	D [м]	U0 [°C]	Θ0 [°C]
0.07	194	0.0022	41	166

 

, где

P [°C]	S [сек-4]	Q [сек-2]
24	1.04	7.93

 

φ1(T)	φ2(Q,T)	f1(Y)	F2(Y)
T4	ln(1+Q·T2)	Y4	(cos(G·Y)·Y)2

 

где        G=0.7

 

Ui[°C]	41.0	53.5	66.1	78.5	97.0	78.5	91.0	103.5	116.0	128.5	166.0
αi[ккал/м2ч°С]	3470	3481	3474	3490	3524	3520	3525	3543	3608	3700	3740

 

 

2. Построение оценки зависимости коэффициента конвективного теплообмена от температуры.

 

Зависимость коэффициента конвективного теплообмена α от температуры пластины U:

U₀ – начальная температура пластины

Θ₀ – температура окружающей среды после момента t_0.

Требуется на основании экспериментальных  данных (U_i, α_i) i=0,1,...,10 построить оценки неизвестных коэффициентов .

Оценки коэффициентов будем искать по методу наименьших квадратов (МНК), т.е. потребуем минимума суммы , где

Необходимое (и в данном случае достаточное) условие минимума S:

 

Эта система линейных алгебраических уравнений с тремя неизвестными ( ) в матричном виде:

P·a = V,  где P = , где P-матрица коэффициентов, a- матрица неизвестных, V- матрица свободных членов.

detP≠0, поэтому система имеет единственное решение a = P^-1·V

 

 

 

 

 

 

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

         XXX- Эмпирическая зависимость коэффициента.

 ___- Зависимость коэффициента по методу наименьших квадратов.

 

Таким образом, оценка зависимости  коэффициента конвективного теплообмена  от температуры будет иметь вид

α(U)=3477.967 +218.287×f₁(U) + 98.015232×f₂(U).

 

 

3. Определение момента установления температуры окружающей среды.

 

Закон изменения окружающей среды:

   где 

Момент времени t₀ определим из уравнения:

где

Согласно исходным данным получаем вид функции:

Отделение корней – поиск, по возможности, тесных промежутков, содержащих единственный корень уравнения. Условия существования  и единственности корня уравнения f(t) = 0 на [a,b]:

1. f(a)·f(b) < 0
2. f(t) – монотонна на [a,b], т.е. или .

Интервал, в котором  находится корень можно определить по графику.

 

Из графика функции  видно, что корень находится в  интервале [1.2;1.4].

 

 

1. Уточнение корня методом половинного деления.

 

Требуется найти приближённое значение корня уравнения с абсолютной погрешностью не превосходящей ε = 0,0001.

 

Для нахождения корня возьмём приближение  , , …

и проверим: если функция f(t) в точках ζ и b принимает значения разных знаков, сдвигаем границу a в точку ζ, а если одного знака, то сдвигаем границу b в точку ζ. Останавливаемся в том случае, если ׀a-b׀<ε.

Результаты вычислений методом  половинного деления:

2. Уточнение корня комбинированным методом (Ньютона).

 

Комбинированный метод  выявляет приближённое значение корня  путём приближения к нему слева  и справа методами хорд и касательных.

 

Из рисунка видно, что метод  хорд мы применяем следующим образом: берём за начальные точки точки A₀(a₀) и B₀(b₀), после чего соединяем их хордой AB. Затем находим точку пересечения этой хорды с осью ОХ и называем её a₁. Эта точка будет первым приближением к корню ζ слева. По методу касательных мы проводим касательную в точке В и также находим точку пересечения касательной с осью ОХ, которую называем b₁. Эта точка будет первым приближением к корню ζ справа.

Полученный отрезок [a₁ b₁] содержится в начальном отрезке [a b], меньше его и тоже содержит корень ζ. Дальнейшие приближения a_i находим тем же способом: проводим хорду между точками A_i-1 и B_i-1. Приближения b_i проводя касательные в точках B_i-1.

Условие применения метода – сохранение выпуклости функции на [a b], т.е. f”(t) > 0 или f “(t) < 0 на всём интервале [a b].

Рассмотрим четыре случая:

 

Из них видно, что  в тех случаях, когда первая и  вторая производная одинаковых знаков (т.е. ) на [a,b], то метод касательных нужно применять насчиная с точки b:

, -метод касательных.

, -метод секущих (хорд).

Если производные разных знаков ( ) на [a,b], то метод касательных(Ньютона) применяется начиная с точки и в формулах мы заменяем на и на . Условие окончания итераций: будем считать, что приближённое значение корня найдено, если вычисления оканчиваются, когда . За корень примем значение:

Результаты, полученные в MathCAD комбинированным методом:

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

3.3 Уточнение корня методом итераций.

Метод итераций применим для решения уравнений специального вида: x = φ(x);

Для решения выбирается начальное приближение x₀, принадлежащее отрезку [a,b]. Последующие приближения находятся следующим образом:

x₁ = φ(x₀), x₂ = φ(x₁), x₃ = φ(x₂), …, x_n = φ(x_n-1)

Доказывается, что последовательность x_n→ζ при n→∞, если ׀φ′(x)׀ ≤ q < 1 на [a,b].

 

 

Условия окончания итераций: , где - требуемая точность. =10^-4.

Уравнение необходимо привести к виду x= φ(x).

 

 

Приведение уравнения f(t) = 0 к виду, пригодному для метода итераций

 

, где 

Назовём то, что в правой части функцией φ(x) и выберем μ  таким образом, чтобы , где В качестве μ возьмём выражение , где m = min f′(t) на отрезке [a,b], а M = max f′(t) на отрезке [a,b].

φ′(t) = 1 - μ·f′(t)

 

1. Допустим m < 0, M < 0

 

2. Допустим m > 0, M > 0

Из пунктов 1 и 2 следует, что 

 

 

Результаты, вычислений в MathCAD методом итераций (формулы и результаты вычислений):

4. Вычисление интеграла I.

 

Для линейной модели с  сосредоточенными параметрами имела  место формула:

,    t ≥ t₀

где

,

Требуется найти приближённое значение интеграла I с абсолютной погрешностью не превосходящей 0,0001(Θ₀ – U₀)

[c]         a₀ – смотри в матрице в п.2

Обозначим через J интеграл I умноженный на T:

Абсолютная погрешность  вычисления интеграла J равна соответственно 0,0001(Θ₀ – U₀)T.

 

1. Вычисление интеграла J по формуле прямоугольников.

 

Дан интеграл , где

За шаг интегрирования выберем величину h = t₀/n, где t₀ – корень f(t) = 0, а n – число разбиений. τ_i = h·i, где i = 0,1,2,…,n;

Точное значение интеграла J = Ĵ + J, где Ĵ – приближённое значение интеграла, а J – погрешность вычисления.

Интеграл считаем как  сумму площадей прямоугольников, т.е.:

Ĵ

Оценка погрешности  формулы прямоугольников:

 Ĵ , где

Зная ε из неравенства    можно получить

Для оценки погрешности  используется правило Рунге:  
если ׀Ĵ_2n – Ĵ_n׀ ≤ 3·ε, тогда ׀J – Ĵ_2n׀ ≤ ε.

 

 

 

Результаты, полученные по формуле прямоугольников с точностью

 T=3.1555.

 

n	Jn	|J2n-Jn|
1		-
2		9.79
4		4.1
8		1.18

 

Значение интеграла  .

2. Вычисление интеграла по формуле трапеций.

Дан интеграл , где

 

Шаг интегрирования выбираем тот же, что и в методе прямоугольников: h = t₀/n;  τ_i = h·i, где i = 0,1,…,n.  Точное значение интеграла J = Ĵ + J

, где Ĵ – приближённое  значение интеграла, вычисленное  по формуле трапеций, а  J – погрешность вычисления.

Приближённое значение интеграла будем считать как  сумму площадей трапеций, показанных на рисунке:

Ĵ =

- формула трапеций

Оценка погрешности  формулы трапеций:

 Ĵ , где

Для оценки J используем правило Рунге:

если  ׀Ĵ_2n – Ĵ_n׀ ≤ 3·ε, тогда ׀J – Ĵ_2n׀ ≤ ε, где ε = 0,0001(Θ₀ – U₀)T.

Результаты, полученные по формуле трапеций с точностью

 T=3.1555.

n	Jn	|J2n-Jn|
1		-
2		27.4986
4		8.85
8		2.37
16		0.6

 

Значение интеграла  .

 

3. Вычисление интеграла J по формуле парабол (Симпсона).

Дан интеграл , где

Число n должно быть обязательно  чётным. Проводим параболы через каждые три рядом стоящие точки. Приближённое значение интеграла считаем по формуле  Симпсона:

Ĵ =

Погрешность формулы  Симпсона:

 Ĵ , где

 

Правило Рунге:

если  ׀Ĵ_2n – Ĵ_n׀ ≤ 15·ε, тогда ׀J – Ĵ_2n׀ ≤ ε, где ε = 0,0001(Θ₀ – U₀)T.

Результаты, полученные по формуле Симпсона с точностью

T=3.1555.

n	Jn	|J2n-Jn|
1		-
2		2.64

Значение интеграла  .

 

5. Приближённое решение задачи Коши.

Для нелинейной модели с  сосредоточенными параметрами было получено уравнение:

- нелинейная з-ча Коши.

,   где 

Сделаем переход к  безразмерной переменным:

Умножим уравнение на :

где

  

 

1. Решение задачи Коши методом Эйлера.

в нашем случае:

Выберем h- шаг интегрирования и построим сетку интегрирования: ; ,   где y_i = y(x_i) – сеточная функция;

;      

Расчётные формулы для  метода Эйлера:

Ломаная Эйлера.

Оценка погрешностей метода Эйлера по правилу Рунге:

, где  y(x_i) – точное значение y(x) в точке x_i; y_i – приближённое значение функции в точке x_i при шаге h; y_2i – приближённое значение функции в точке x_i при шаге h/2.

Результаты, полученные в MathCad’е:

T=3.1555 b1=6.2763·10-2 b2=2.8182·10-2

P=0.192 X0=0.3986

i	xi	yi	y2i	|y2i- yi|
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15	0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0	0 0 6.1061·10-2 0.2489 0.3994 0.5202 0.6172 0.6952 0.7577 0.8077 0.8477 0.8796 0.9049 0.925 0.9409 0.9535	0 1.1372·10-2 9.2908·10-2 0.2654 0.4053 0.5191 0.6116 0.6867 0.7477 0.7971 0.8371 0.8671 0.8954 0.9164 0.9332 0.9466	0 1.1372 ·10-2 0.031847 0.0165 0.0059 0.0011 0.0056 0.0085 0.01 0.0106 0.0106 0.0125 0.0095 0.0086 0.0077 0.0069

 

2. Решение задачи Коши методом Рунге-Кута.

в нашем случае:

Расчётные формулы метода Рунге-Кута.

Пусть h – шаг интегрирования.

где  

i = 0,1,2,...

Оценка погрешности  метода Рунге-Кута:

y(x_i) – точное значение y(x) в точке x_i;

y_2i – приближённое значение y(x_i), полученное при шаге интегрирования h/2;

y_i – приближённое значение y(x_i), полученное при шаге интегрирования h/2;

0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0	0 2.39·10-2 0.1291 0.2872 0.4169 0.5234 0.6109 0.6828 0.7417 0.7899 0.8294 0.8616 0.8878 0.9092 0.9265 0.9406	0 2.3581·10-2 0.1288 0.287 0.4167 0.5233 0.6108 0.6827 0.7416 0.7899 0.8294 0.8616 0.8878 0.9092 0.9265 0.9406	0 2.126·10-5 2·10-5 1.44·10-4 1.33·10-5 6.66·10-6 6.66·10-6 6.66·10-6 6.66·10-6 0 0 0 0 0 0 0

 

 

Возвращение к размерным  переменным.

Для того, чтобы построить  интересующую нас зависимость U(t) для нелинейной модели, необходимо вернуться к размерным переменным по формулам:

t_i = x_i∙T [c],   U_i = U₀ + (Θ₀ - U₀)∙y_I [ºC],

 где данные взяты по методу Рунге-Кутта при шаге h=0.2.

Для линейной модели с  сосредоточенными параметрами:

,

где интеграл I найден при помощи метода парабол.

ti	Unelini	Ulini
0	41	-
0.6311	43.9875	-
1.2622	57.1349	57.1652
1.8933	76.8941	76.8936
2.5244	93.1102	93.0458
3.1555	106.35	106.277
3.7866	117.37	117.1
4.4178	126.356	125.96
5.0489	133.712	133.22
5.68	139.745	139.16
6.3111	144.689	144.03
6.9422	148.7	148.01
7.5733	151.986	151.27
8.2044	154.65	153.94
8.8355	156.82	156.13
9.4666	158.58	157.92