Санкт-Петербургский Государственный Технологический Институт
(Технический Университет)
Кафедра Факультет VIII
Прикладной Курс II
Математики Группа 891
Дисциплина: Информатика – 2
Курсовая работа
Тема: “Исследование распределения температуры в тонком цилиндрическом стержне”
Руководитель:
Поляков В.О.
Исполнитель:
Солнцев П.В.
Санкт-Петербург 2001
Введение
В решении любой прикладной задачи можно выделить три основных этапа:
-
построение математической модели исследуемого объекта
-
выбор способа и алгоритма решения полученной модели
-
численная реализация алгоритма
Цель данной работы – на примере исследования распределения температуры в тонком цилиндрическом стержне освоить основные методы приближённых вычислений, приобрести практические навыки самостоятельных исследований, существенно опирающихся на использование методов прикладной математики.
Содержание
-
Постановка задачи
-
Физическая модель
-
Математическая модель
-
Обработка результатов эксперимента
-
Задача регрессии. Метод наименьших квадратов.
-
Гипотеза об адекватности модели задачи регрессии
-
Нахождение коэффициента теплоотдачи a
-
Вычисление интеграла методом трапеций
-
Вычисление интеграла методом парабол (Симпсона)
-
Вычисление времени Т0 установления режима
-
Решение уравнения комбинированным методом
-
Решение уравнения методом итерраций
-
Решение краевой задачи (метод малого параметра)
-
Заключение
Литература
-
Постановка задачи
-
Физическая модель
В ряде практических задач возникает необходимость исследования распределения температуры вдоль тонкого цилиндрического стержня, помещённого в высокотемпературный поток жидкости или газа. Это исследование может проводиться либо на основе обработки эксперимента (измерение температуры в различных точках стержня), либо путём анализа соответствующей математической модели.
В настоящей работе используются оба подхода.
Тонкий цилиндрический стержень помещён в тепловой поток с постоянной температурой q, на концах стержня поддерживается постоянная температура q0.
1.2 Математическая модель
Совместим координатную ось абсцисс с продольной осью стержня с началом в середине стержня. Будем рассматривать задачу (распределения температуры по стержню) мосле момента установления режима Т0.
Первая математическая модель использует экспериментальные данные, при этом измеряют температуру Ui стержня в нескольких точках стержня с координатами xi. Результаты измерения Ui рассматривают как функцию регрессии и получают статистики. Учитывая чётность U(x) можно искать её в виде многочлена по чётным степеням x (ограничимся 4-ой степенью этого многочлена).
(1.1)
Задача сводится к отысканию оценок неизвестных параметров, т.е. коэффициентов a0 , a1 и a2 , например, методом наименьших квадратов.
Вторая математическая модель, также использующая экспериментальные данные, состоит в применении интерполяционных формул и может употребляться, если погрешность измерений температуры Ui пренебрежимо мала, т.е. можно считать, что U(xi)=Ui
Третья математическая модель основана на использовании закона теплофизики. Можно доказать, что искомая функция U(x) имеет вид:
(1.2)
где l - коэффициент теплопроводности, a - коэффициент теплоотдачи, D – диаметр стержня, q - температура потока, в который помещён стержень.
Ищем U(x) как решение краевой задачи для уравнения (1.2) с граничными условиями:
(1.3)
на отрезке [-L|/2;L/2], где L – длина стержня, q0 - постоянная температура, поддерживаемая на концах стержня.
Коэффициент теплопроводности l зависит от температуры:
(1.4)
где l0 - начальное значение коэффициента теплопроводности, sl - вспомогательный коэффициент.
Коэффициент теплоотдачи a вычисляют по формуле:
(1.5)
т.е. как среднее значение функции
за некоторый отрезок времени от 0 до Т, здесь a0 - значение a при t стремящемся к бесконечности, b – известный коэффициент.
Время Т0, по истечении которого распределение температуры в стержне можно считать установившимся определяется по формуле:
(1.6)
где а – коэффициент температуропроводности, x - наименьший положительный корень уравнения:
(1.7)
Задание курсовой работы
Вариант № 136
Исходные данные:
-
L = 0.0386 м
-
D = 0,00386 м
-
q = 740 оС
-
q0 = 74 оС
-
l0 = 141,85 (Вт/м*К)
-
sl = 2,703*10-4
-
B = 6,789*10-7
-
a0 = 3,383*102 (Вт/м2*К)
-
T = 218 оС
-
А = 3,043*10-5 (м2/с)
11
X, м | U, oC |
0 | 353 |
0,00386 | 343 |
0,00772 | 313 |
0,01158 | 261 |
0,01544 | 184 |
0,01930 | 74 |
2. Обработка результатов эксперимента.
2.1 Задача регрессии. Метод наименьших квадратов.
Ищем функцию регрессии в виде (1.1). Оценки коэффициентов находим с помощью МНК, при этом наименьшими будут оценки, обеспечивающие минимум квадратов отклонений оценочной функции регрессии от экспериментальных значений температуры; суммирование ведут по всем экспериментальным точкам, т.е. минимум величины S:
(2.1)
В нашем случае необходимым т достаточным условием минимума S будут:
Где k = 0, 1, 2. (2,2)
Из уравнений (2.1) и (2.2) получаем:
(2.3)
Сумма
Система (2.3) примет вид:
(2.4)
В результате вычислений получаем Sk и Vj. Обозначим матрицу коэффициентов уравнения (2.4) через “p”:
Методом Гаусса решаем систему (2.4) и найдём обратную матрицу p-1. В результате получаем:
Подставляя в (2.1) найденные значения оценок коэффициентов ак, находим минимальное значение суммы S:
Smin=0.7597
При построении доверительных интервалов для оценок коэффициентов определяем предварительно точечные оценки.
Предполагается, что экспериментальные значения xi измерены с пренебрежимо малыми ошибками, а случайные ошибки измерения величины Ui независимы и распределены по нормальному закону с постоянной дисперсией s2, которая неизвестна. Для имеющихся измерений температуры Ui неизвестная дисперсия оценивается по формуле:
Где r – число степеней свободы системы, равное разности между количеством экспериментальных точек и количеством вычисляемых оценок коэффициентов, т.е. r = 3.
Оценка корреляционной матрицы имеет вид:
Оценки дисперсий параметров оценок коэффициентов найдём по формулам:
Где Sk – минор соответствующего диагонального элемента матрицы нормальной системы;
D - главный определитель нормальной системы.
В нашем случае:
S0=3.5438 10-22
S1=-8.9667 10-14
S2=6.3247 10-7
Откуда:
Найденные оценки коэффициентов распределены по нормальному закону, т.к. линейно зависят от линейно распределённых экспериментальных данных Ui.
Известно, что эти оценки несмещённые и эффективные. Тогда случайные величины:
Имеют распределения Стьюдента, а r = 3.
Выбираем доверительную вероятность b=0,9 и по таблице Стьюдента находим критическое значение gb равное 2,35, удовлетворяющее равенству:
Доверительные интервалы для коэффициентов:
(2.4*)
В нашем случае примут вид:
2.2 Проверка статистической гипотезы об адекватности модели задачи регрессии.
Имеется выборка объёма n экспериментальных значений (xi;Ui). Предполагаем, что ошибки измерения xi пренебрежимо малы, а случайные ошибки измерения температур Ui подчинены нормальному закону с постоянной дисперсией s2. Мы выбрали функцию регрессии в виде:
Выясним, нельзя ли было ограничиться многочленом второго порядка, т.е. функцией вида:
(2.5)
C помощью МНК можно найти оценки этих функций и несмещённый оценки дисперсии отдельного измерения Ui для этих случаев:
Где r1 = 4 (количество точек – 6, параметра – 2).
Нормальная система уравнений для определения новых оценок коэффициентов функции (2.5)с помощью МНК имеет вид:
(2.7)
Решая эту систему методом Гаусса, получим:
(2.8)
Чем лучше функция регрессии описывает эксперимент, тем меньше для неё должна быть оценка дисперсии отдельного измерения Ui, т.к. при плохом выборе функции в дисперсию войдут связанные с этим выбором дополнительные погрешности. Поэтому для того, чтобы сделать выбор между функциями U(x) и U(1)(x) нужно проверить значимость различия между соответствующими оценками дисперсии, т.е. проверить гипотезу:
Н0 – альтернативная гипотеза
Т.е. проверить, значимо ли уменьшение дисперсии при увеличении степени многочлена.
В качестве статического критерия рассмотрим случайную величину, равную:
(2.9)
имеющую распределение Фишера с(r ; r1) степенями свободы. Выбираем уровень распределения Фишера, находим критическое значение F*a, удовлетворяющее равенству: p(F>F*a)=a
В нашем случае F=349.02, а F*a=10,13.
Если бы выполнилось практически невозможное соотношение F>Fa, имевшее вероятность 0,01, то гипотезу Н0 пришлось бы отклонить. Но в нашем случае можно ограничиться многочленом
, коэффициенты в котором неодинаковы.
3. Нахождение коэффициента теплопроводности a.
Коэффициент a вычислим по формуле (1.5), обозначим:
(3.1)
Определим допустимую абсолютную погрешность величины интеграла I, исходя из требования, чтобы относительная погрешность вычисления a не превосходила 0,1%, т.е.:
(3.2)
Т.к. из (3.1) очевидно, что a>a0, то условие (3.2) заведомо будет выполнено, если:
(3.3)
Т.е. в качестве предельно допустимой абсолютной погрешности вычисления интеграла I возьмём d=0,001Т (3.4)
Т=218 оС, следовательно, d=0,218 оС.
3.1 Вычисление интеграла I методом трапеции
Использование теоретической оценки погрешности
Для обозначения требуемой точности количества частей n, на которые нужно разбить отрезок интегрирования [0;T] определяется по формуле:
, где M2=[f”(t)], t e [0;T], f(t)=e-bt3
Учитывая формулу (3.4) получаем:
(3.5)
Дифференцируя f(t), получим:
А необходимое условие экстремума: f”(t)-f’’’(t)=0, откуда получаем:
Далее вычисляем значения f’’(t) при t=t1, t=t2, t=0 и t=T, получаем:
f’’(t1)=1.5886 10-4
f’’(t2)=-1.6627 10-4
f’’(0)=0
f’’(T)=7.4782 10-6
Итак: M2=1,5886 10-4, откуда n=25.66; принимаем N=26.
Далее вычислим интеграл I:
Погрешность вычисления a:
3.2 Вычисление интеграла I методом парабол
При расчётах будем использовать теоретическую оценку погрешности с помощью правила Рунге. Для обеспечения заданной точности количество частей n, на которое следует разделить интервал интегрирования можно определить по формуле:
, откуда:
Нахождение М4 можно провести аналогично нахождению М2 в предыдущем пункте, но выражение для fIV(t) имеет довольно громоздкий вид. Поэтому правило Рунге – наиболее простой способ.
Обозначим через In и I2n значение интеграла I, полученное при разбиении промежутка интегрирования соответственно на n и 2n интервалов. Если выполнено равенство: |I2n-In| = 15d (*1), то |I-I2n|=d
Будем , начиная с n=2, удваивать n до тех пор, пока не начнёт выполняться неравенство (*1), тогда:
(3.6)
Согласно формуле парабол (3.7):
Результаты вычислений сведём в таблицу:
n | In | I2n |
4 | 102.11 | |
8 | 101.61 | 0.5017 |
По формуле (3.7) I = 101,61 что в пределах погрешности совпадает со значением, полученным по методу трапеций
n=8 | n=4 |
ti (8) | y8 | ti (4) | y4 |
0 | 1 | 0 | 1 |
27.25 | 0.9864 | | |
54.5 | 0.8959 | 54.5 | 0.8959 |
81.75 | 0.6901 | | |
109 | 0.4151 | 109 | 0.4151 |
136.25 | 0.1796 | | |
163.5 | 0.0514 | 163.5 | 0.0514 |
190.75 | 0.0089874 | | |
218 | 0.00088179 | 218 | 0.00088179 |
4. Вычисление времени Т0 установления режима
4.1 Решение уравнения комбинированным методом
Время установления режима определяется по формулам (1.6) и (1.7).
Проведём сначала отделение корней. Имеем y = ctg(x) и y = Ax. Приведём уравнение к виду: A x sin(x)-cos(x) = 0. Проведём процесс отделения корня.
F(x) | -1 | -0.6285 | 0.4843 |
x | 0.01 | 0.05 | 0.1 |
т.е. x с [0.01;0.05]
Убедимся, что корень действительно существует и является единственным на выбранном интервале изоляции.
f(a) f(b)<0 – условие существования корня выполняется
f’(x) на [a;b] – знакопостоянна: f’(x)>0 – условие единственности также выполняется. Проведём уточнение с погрешностью не превышающей e=10-4
Строим касательные с того конца, где f(x) f”(x)>0
f”(x)=(2A+1)cos(x) – A x sin(x). f”(x)>0 на (a;b), следовательно касательные строим справа, а хорды слева. Приближение корня по методу касательных:
по методу хорд:
Вычисление ведём до того момента, пока не выполнится условие:
Результаты вычислений заносим в таблицу:
n | an | bn | f(an) | f(bn) |
0 | 0.05 | 0.1 | -0.6285 | 0.4843 |
1 | 0.07824 | 0.08366 | -0.0908 | 0.0394 |
2 | 0.08202 | 0.08207 | -9.1515 10-4 | 3.7121 10-4 |
3 | 0.08206 | 0.08206 | -8.4666 10-8 | 3.4321 10-8 |
Т0 = 72,7176 секунд.
4.2 Решение уравнения комбинированным методом
Приведём f(x) = 0 к виду x = j(x). Для этого умножим обе части на произвольное число m, неравное нулю, и добавим к обеим частям х:
X = x - m f(x)
j(x) = x - m A x sin(x) + m cos(x)
В качестве m возьмём:
где М = max [f’(x)] на [a;b], а m = min [f’(x)] на [a’b]
В силу монотонности f’(x) на [a;b] имеем m = f’(а), М = f’(b). Тогда m = 0,045.
Приближение к корню ищем по следующей схеме:
Вычисление ведём до тех пор, пока не выполнится условие:
(q = max |j’(x)| на [a’b])
j’(x) на [a’b] монотонно убывает, поэтому максимум его модуля достигается на одном из концов.
j’(0,05) = 0,3322 j’(0,1) = -0,3322, следовательно, q = 0.3322 < 1. В этом случае выполняется условие сходимости и получается последовательность:
i | xi | j( xi) | D xi |
0 | 0.075 | 0.082392 | 0.00739 |
1 | 0.082392 | 0.082025 | 0.000367 |
2 | 0.082025 | 0.08206 | 3.54 10-5 |
3 | 0.08206 | 0.082057 | 3.33 10-6 |
4 | 0.082057 | 0.082057 | 3.15 10-7 |
Итак, с погрешностью, меньшей 10-4, имеем:
Т0 = 72,7176 с. , x = 0.03142
5. Решение краевой задачи
Используем метод малого параметра. Краевую задачу запишем в виде:
(5.1)
Введя новую переменную y = (U - q0)/(q - q0), запишем (5.1) в виде:
(5.2)
e = sl(q - q0) =0.18, L/2 =0.0193. В качестве малого параметра возьмём e.
Тогда, подставив y(x) в уравнение (5.2) и перегруппировав члены при одинаковых степенях e, получим:
(5.3)
Ограничимся двумя первыми членами ряда:
Из (5.2) и (5.3) находим общее решение уравнения для y0:
где y0 с тильдой – частное решение данного неоднородного уравнения; y(1) и y(2) – линейно независимые решения однородного уравнения.
Корни уравнения:
y0общ = 1 + c1ch(px)+c2sh(px), где p = 0.01953
Константы найдём из граничных условий:
откуда с1 = 0, с2 = -0,57; т.е. имеем функцию:
y0 = 1 - 0.57 sh(px)
Общее решение:
Частное решение:
Дифференцируя и подставляя в уравнение, получим:
А1 = 0; А2 = -0,1083; В1 = 0; В2 = 17,1569;
Тогда общее решение для y1 имеет вид:
с3 = 0; с4 = 0,0462
Перейдя к старой переменной U, получим:
q0 = 0; q1 = -374.11; q2 = -12.9863; q3 = 2057
Итоговое уравнение:
Пользуясь этой формулой, составим таблицу значений функции U(x):
x | U(x) | U |
0 | 352.9075 | 353 |
0.0019 | 350.4901 | |
0.0039 | 343.1972 | 343 |
0.0058 | 330.9053 | |
0.0077 | 313.4042 | 313 |
0.0097 | 290.391 | |
0.0116 | 261.4598 | 261 |
0.0135 | 226.0893 | |
0.0154 | 1836255 | 184 |
0.0174 | 133.2579 | |
0.0193 | 74 | 74 |
Используя данную таблицу, строим график функции U(x).
[см. приложение 1]
6. Заключение
Решение задачи на ЭВМ при помощи вычислительной системы ManhCad 7.0 дало результаты (функцию распределения температуры в тонком цилиндрическом стержне), полученные по решению практического задания и обработкой эксперимента (функции регрессии), которые практически (в пределах погрешности) совпадают с экспериментальными значениями.
Литература
1. Методические указания “Методы приближённых вычислений. Решение нелинейных уравнений”
(ЛТИ им. Ленсовета, Л. 1983)
2.Методические указания “Приближённые методы ислисления определённых интегралов”
(ЛТИ им. Ленсовета, Л. 1986)
-
Методические указания “Изучение распределения температуры в тонком цилиндрическом стержне”
(ЛТИ им. Ленсовета, Л. 1988)
Приложение 1