Файл: Расчет нестационарного температурного поля в двухслойной бесконечной пластине (изолированном стержне) с заданными начальными и граничными условиями.docx
ВУЗ: Не указан
Категория: Не указан
Дисциплина: Не указана
Добавлен: 27.04.2024
Просмотров: 9
Скачиваний: 0
ВНИМАНИЕ! Если данный файл нарушает Ваши авторские права, то обязательно сообщите нам.
и 
необходимо знать 
и 
которые находятся из левого граничного условия.
Далее по формулам (6) последовательно находятся
, при условии, что 
найдено из правого граничного условия.
Таким образом, решение уравнений вида (5) описываемым способом, называемым методом прогонки, сводится к вычислениям по трем формулам: нахождение так называемых прогоночных коэффициентов
и по формулам (7) i = 2,3,...,N-1 (прямая прогонка) и затем получение неизвестных 
по формуле (6) при i = N-1, N-2,...,2 (обратная прогонка).
Для успешного применения метода прогонки нужно, чтобы в процессе вычислений не возникло ситуаций с делением на нуль, а при больших размерностях систем не должно быть быстрого роста погрешностей округлений.
Будем называть прогонку корректной, если знаменатели прогоночных коэффициентов (7) не обращаются в нуль, и устойчивой, если
при всех 
.
Возвращаясь к системе (4), определим прогоночные коэффициенты и воссоздадим полный алгоритм решения полученной системы.
Поскольку при х = 0 Т = Тл,
то
,
, 
а при
х = L Т = Тn
.
Прогоночные коэффициенты вычисляются по формулам (7).
2.2 Конечно разностная аппроксимация граничных условий
Будем использовать неявную разностную схему, а граничные условия будут использоваться для нахождения первых прогоночных коэффициентов α
1 и β1 (левое граничное условие), и температуры ТN (правое граничное условие).
Граничное условие – 1-го рода:
,
.
Граничное условие – 2-го рода:
.
Граничное условие – 3-го рода:
Граничное условие – 4-го рода, определяет прогоночные коэффициенты на границе раздела двух материалов:
3. Основная часть
3.1 Исходные данные для расчета
Рассмотрим краевую задачу на основе одномерного уравнения теплопроводности. Анализируется теплопередача через двухслоную бесконечную пластину (или изолированный стержень). Геометрия задачи показана на рис. 1. На одной стороне пластины (слева) поддерживается температураT=70°С. На другой стороне пластины (справа) задана q=0.
При заданных условиях температура будет изменяться только в направлении оси Ox, теплофизические коэффициенты не зависят от температуры.
3.2 Математическая постановка задачи
Запишем дифференциальные уравнения описывающие процесс теплопередачи в пластине:
Начальные и граничные условия запишем следующим образом
.
Код программы :
Private Sub CommandButton1_Click()
Const mf = 1000
Const pi = 3.141592
Dim T(1 To mf) As Double
Dim Tt(1 To mf) As Double
Dim alfa(1 To mf) As Double
Dim beta(1 To mf) As Double
Dim i, j, N As Double
Dim ai, bi, ci, fi As Double
Dim lamda1, ro1, c1, a1, h, tau, d, Isv, Sd, p_ud1, p_ud2, Q_contact As Double
Dim lamda2, ro2, c2, a2 As Double
Dim T0, Tr, L1, L2, t_end, time As Double
Dim left, right, Ng As Integer
Dim T1, q1, Te1, k1, Qintern1, Qintern As Double
Dim T2, q2, Te2, k2, Qintern2 As Double
N = CSng(TextBox1.Text)
t_kon = CSng(TextBox2.Text)
t_svar = CSng(TextBox29.Text)
L1 = CSng(TextBox21.Text)
L2 = CSng(TextBox22.Text)
lamda1 = CSng(TextBox4.Text)
ro1 = CSng(TextBox5.Text)
c1 = CSng(TextBox6.Text)
lamda2 = CSng(TextBox20.Text)
ro2 = CSng(TextBox18.Text)
c2 = CSng(TextBox19.Text)
T0 = CSng(TextBox7.Text)
T1 = CSng(TextBox10.Text)
T2 = CSng(TextBox12.Text)
Te1 = CSng(TextBox15.Text)
Te2 = CSng(TextBox14.Text)
q1 = CSng(TextBox11.Text)
q2 = CSng(TextBox13.Text)
k1 = CSng(TextBox16.Text)
k2 = CSng(TextBox17.Text)
d = CSng(TextBox25.Text)
Isv = CSng(TextBox24.Text)
a1 = lamda1 / (c1 * ro1)
a2 = lamda2 / (c2 * ro2)
Sd = pi * d ^ 2 / 4
j = Isv / Sd
p_ud1 = CSng(TextBox27.Text)
p_ud2 = CSng(TextBox26.Text)
Qintern1 = p_ud1 * j ^ 2
Qintern2 = p_ud2 * j ^ 2
Q_contact = CSng(TextBox28.Text)
'шаг сетки по пространственной координате
h = (L1 + L2) / (N - 1)
Ng = L1 / h
'расчетный шаг сетки по времени
tau = t_kon / 100
'Определяем поле температур в начальный момент времени
For i = 1 To N
T(i) = T0
Next i
If OptionButton1.Value = True Then
left = 1
End If
If OptionButton2.Value = True Then
left = 2
End If
If OptionButton3.Value = True Then
left = 3
End If
If OptionButton4.Value = True Then
right = 1
End If
If OptionButton5.Value = True Then
right = 2
End If
If OptionButton6.Value = True Then
right = 3
End If
'Интегрируем нестационарное уравнение теплопроводности
time = 0
Dim time_i As Integer
time_i = 0
While time < t_kon
time = time + tau
time_i = time_i + 1
If time > t_svar Then
Q_contact = 0
Qintern = 0
Qintern1 = 0
Qintern2 = 0
End If
'Определяем начальные прогоночные коэффициенты на основе левого граничного условия
Select Case left
Case 1
alfa(1) = 0
beta(1) = T1
Case 2
alfa(1) = (2 * a1 * tau) / (h ^ 2 + 2 * a1 * tau)
beta(1) = (h ^ 2 / (h ^ 2 + 2 * a1 * tau)) * T(1) + (2 * a1 * tau * h * q1) / lamda1 * (h ^ 2 + 2 * a1 * tau)
Case 3
alfa(1) = (2 * a1 * tau * lamda1) / (lamda1 * h ^ 2 + 2 * a1 * tau * (lamda1 + h * k1))
beta(1) = (lamda1 * h ^ 2 * T(1) + 2 * a1 * tau * h * k1 * Te1) / (lamda1 * h ^ 2 + 2 * a1 * tau * (lamda1 + h * k1))
End Select
'Определяем прогоночные коэффициенты
For i = 2 To Ng - 1
If i = Ng - 1 Then
Qintern = Q_contact
Else
Qintern = Qintern1
End If
ai = lamda1 / h ^ 2
bi = 2 * lamda1 / h ^ 2 + ro1 * c1 / tau
ci = lamda1 / h ^ 2
fi = -ro1 * c1 * T(i) / tau - Qintern
alfa(i) = ai / (bi - ci * alfa(i - 1))
beta(i) = (ci * beta(i - 1) - fi) / (bi - ci * alfa(i - 1))
Next i
alfa(Ng) = 2 * a1 * a2 * tau * lamda2 / (2 * a1 * a2 * tau * (lamda2 + lamda1 * (1 - alfa(Ng - 1))) + h ^ 2 * (-a1 * lamda1 + a2 * lamda2))
beta(Ng) = (2 * a1 * a2 * tau * lamda1 * beta(Ng - 1) + h ^ 2 * (-a1 * lamda1 + a2 * lamda2) * T(Ng)) / (2 * a1 * a2 * tau * (lamda2 + lamda1 * (1 - alfa(Ng - 1))) + h ^ 2 * (-a1 * lamda1 + a2 * lamda2))
For i = Ng + 1 To N - 1
ai = lamda2 / h ^ 2
bi = 2 * lamda2 / h ^ 2 + ro2 * c2 / tau
ci = lamda2 / h ^ 2
fi = -ro2 * c2 * T(i) / tau - Qintern2
alfa(i) = ai / (bi - ci * alfa(i - 1))
beta(i) = (ci * beta(i - 1) - fi) / (bi - ci * alfa(i - 1))
Next i
'Определяем значение температуры на правой границе
Select Case right
Case 1
T(N) = T2
Case 2
T(N) = (2 * a2 * tau * lamda2 * beta(N - 1) - 2 * a2 * tau * h * q2 + h ^ 2 * lamda2 * T(N)) / (lamda2 * h ^ 2 + 2 * a2 * tau * lamda2 * (1 - alfa(N - 1)))
Case 3
T(N) = (lamda2 * h ^ 2 * T(N) + 2 * a2 * tau * (lamda2 * beta(N - 1) + h * k2 * Te2)) / (lamda2 * h ^ 2 + 2 * a2 * tau * (h * k2 + lamda2 * (1 - alfa(N - 1))))
End Select
'Определяем поле температур
For i = (N - 1) To 1 Step -1
T(i) = alfa(i) * T(i + 1) + beta(i)
Next i
Tt(time_i) = T(Ng - 1)
Wend
'Записываем результаты T(i) в таблицу
Dim S As String 'пользовательская переменная
Dim v_Excel As Excel.Application 'область памяти приложения Excel
Dim v_Wb1 As Excel.Workbook 'область рабочей книги
Dim sh As Excel.Worksheet 'область рабочего листа
Set v_Excel = New Excel.Application 'Открываем новое приложение Excel
'Загружаем рабочую книгу
Set v_Wb1 = v_Excel.Workbooks.Open("I:temp.xlsx")
Set sh = v_Excel.Sheets(1) 'Загружаем рабочий лист 1
v_Wb1.Activate 'Получаем доступ к рабочей книге
For i = 1 To N
sh.Cells(i, 1).Value = (i - 1) * h
sh.Cells(i, 2).Value = T(i)
Next i
For i = 1 To 100
sh.Cells(i, 3).Value = (i - 1) * tau
sh.Cells(i, 4).Value = Tt(i)
Next i
v_Excel.Quit 'Закрытие приложения Excel
Set v_Excel = Nothing
End Sub
необходимо знать и которые находятся из левого граничного условия.Далее по формулам (6) последовательно находятся
, при условии, что найдено из правого граничного условия.Таким образом, решение уравнений вида (5) описываемым способом, называемым методом прогонки, сводится к вычислениям по трем формулам: нахождение так называемых прогоночных коэффициентов
и по формулам (7) i = 2,3,...,N-1 (прямая прогонка) и затем получение неизвестных по формуле (6) при i = N-1, N-2,...,2 (обратная прогонка).Для успешного применения метода прогонки нужно, чтобы в процессе вычислений не возникло ситуаций с делением на нуль, а при больших размерностях систем не должно быть быстрого роста погрешностей округлений.
Будем называть прогонку корректной, если знаменатели прогоночных коэффициентов (7) не обращаются в нуль, и устойчивой, если
при всех .Возвращаясь к системе (4), определим прогоночные коэффициенты и воссоздадим полный алгоритм решения полученной системы.
Поскольку при х = 0 Т = Тл,
то
, , а при
х = L Т = Тn
.Прогоночные коэффициенты вычисляются по формулам (7).
2.2 Конечно разностная аппроксимация граничных условий
Будем использовать неявную разностную схему, а граничные условия будут использоваться для нахождения первых прогоночных коэффициентов α
1 и β1 (левое граничное условие), и температуры ТN (правое граничное условие).
Граничное условие – 1-го рода:
, .Граничное условие – 2-го рода:
.Граничное условие – 3-го рода:
Граничное условие – 4-го рода, определяет прогоночные коэффициенты на границе раздела двух материалов:
3. Основная часть
3.1 Исходные данные для расчета
Рассмотрим краевую задачу на основе одномерного уравнения теплопроводности. Анализируется теплопередача через двухслоную бесконечную пластину (или изолированный стержень). Геометрия задачи показана на рис. 1. На одной стороне пластины (слева) поддерживается температураT=70°С. На другой стороне пластины (справа) задана q=0.
При заданных условиях температура будет изменяться только в направлении оси Ox, теплофизические коэффициенты не зависят от температуры.
3.2 Математическая постановка задачи
Запишем дифференциальные уравнения описывающие процесс теплопередачи в пластине:
Начальные и граничные условия запишем следующим образом
. Код программы :
Private Sub CommandButton1_Click()
Const mf = 1000
Const pi = 3.141592
Dim T(1 To mf) As Double
Dim Tt(1 To mf) As Double
Dim alfa(1 To mf) As Double
Dim beta(1 To mf) As Double
Dim i, j, N As Double
Dim ai, bi, ci, fi As Double
Dim lamda1, ro1, c1, a1, h, tau, d, Isv, Sd, p_ud1, p_ud2, Q_contact As Double
Dim lamda2, ro2, c2, a2 As Double
Dim T0, Tr, L1, L2, t_end, time As Double
Dim left, right, Ng As Integer
Dim T1, q1, Te1, k1, Qintern1, Qintern As Double
Dim T2, q2, Te2, k2, Qintern2 As Double
N = CSng(TextBox1.Text)
t_kon = CSng(TextBox2.Text)
t_svar = CSng(TextBox29.Text)
L1 = CSng(TextBox21.Text)
L2 = CSng(TextBox22.Text)
lamda1 = CSng(TextBox4.Text)
ro1 = CSng(TextBox5.Text)
c1 = CSng(TextBox6.Text)
lamda2 = CSng(TextBox20.Text)
ro2 = CSng(TextBox18.Text)
c2 = CSng(TextBox19.Text)
T0 = CSng(TextBox7.Text)
T1 = CSng(TextBox10.Text)
T2 = CSng(TextBox12.Text)
Te1 = CSng(TextBox15.Text)
Te2 = CSng(TextBox14.Text)
q1 = CSng(TextBox11.Text)
q2 = CSng(TextBox13.Text)
k1 = CSng(TextBox16.Text)
k2 = CSng(TextBox17.Text)
d = CSng(TextBox25.Text)
Isv = CSng(TextBox24.Text)
a1 = lamda1 / (c1 * ro1)
a2 = lamda2 / (c2 * ro2)
Sd = pi * d ^ 2 / 4
j = Isv / Sd
p_ud1 = CSng(TextBox27.Text)
p_ud2 = CSng(TextBox26.Text)
Qintern1 = p_ud1 * j ^ 2
Qintern2 = p_ud2 * j ^ 2
Q_contact = CSng(TextBox28.Text)
'шаг сетки по пространственной координате
h = (L1 + L2) / (N - 1)
Ng = L1 / h
'расчетный шаг сетки по времени
tau = t_kon / 100
'Определяем поле температур в начальный момент времени
For i = 1 To N
T(i) = T0
Next i
If OptionButton1.Value = True Then
left = 1
End If
If OptionButton2.Value = True Then
left = 2
End If
If OptionButton3.Value = True Then
left = 3
End If
If OptionButton4.Value = True Then
right = 1
End If
If OptionButton5.Value = True Then
right = 2
End If
If OptionButton6.Value = True Then
right = 3
End If
'Интегрируем нестационарное уравнение теплопроводности
time = 0
Dim time_i As Integer
time_i = 0
While time < t_kon
time = time + tau
time_i = time_i + 1
If time > t_svar Then
Q_contact = 0
Qintern = 0
Qintern1 = 0
Qintern2 = 0
End If
'Определяем начальные прогоночные коэффициенты на основе левого граничного условия
Select Case left
Case 1
alfa(1) = 0
beta(1) = T1
Case 2
alfa(1) = (2 * a1 * tau) / (h ^ 2 + 2 * a1 * tau)
beta(1) = (h ^ 2 / (h ^ 2 + 2 * a1 * tau)) * T(1) + (2 * a1 * tau * h * q1) / lamda1 * (h ^ 2 + 2 * a1 * tau)
Case 3
alfa(1) = (2 * a1 * tau * lamda1) / (lamda1 * h ^ 2 + 2 * a1 * tau * (lamda1 + h * k1))
beta(1) = (lamda1 * h ^ 2 * T(1) + 2 * a1 * tau * h * k1 * Te1) / (lamda1 * h ^ 2 + 2 * a1 * tau * (lamda1 + h * k1))
End Select
'Определяем прогоночные коэффициенты
For i = 2 To Ng - 1
If i = Ng - 1 Then
Qintern = Q_contact
Else
Qintern = Qintern1
End If
ai = lamda1 / h ^ 2
bi = 2 * lamda1 / h ^ 2 + ro1 * c1 / tau
ci = lamda1 / h ^ 2
fi = -ro1 * c1 * T(i) / tau - Qintern
alfa(i) = ai / (bi - ci * alfa(i - 1))
beta(i) = (ci * beta(i - 1) - fi) / (bi - ci * alfa(i - 1))
Next i
alfa(Ng) = 2 * a1 * a2 * tau * lamda2 / (2 * a1 * a2 * tau * (lamda2 + lamda1 * (1 - alfa(Ng - 1))) + h ^ 2 * (-a1 * lamda1 + a2 * lamda2))
beta(Ng) = (2 * a1 * a2 * tau * lamda1 * beta(Ng - 1) + h ^ 2 * (-a1 * lamda1 + a2 * lamda2) * T(Ng)) / (2 * a1 * a2 * tau * (lamda2 + lamda1 * (1 - alfa(Ng - 1))) + h ^ 2 * (-a1 * lamda1 + a2 * lamda2))
For i = Ng + 1 To N - 1
ai = lamda2 / h ^ 2
bi = 2 * lamda2 / h ^ 2 + ro2 * c2 / tau
ci = lamda2 / h ^ 2
fi = -ro2 * c2 * T(i) / tau - Qintern2
alfa(i) = ai / (bi - ci * alfa(i - 1))
beta(i) = (ci * beta(i - 1) - fi) / (bi - ci * alfa(i - 1))
Next i
'Определяем значение температуры на правой границе
Select Case right
Case 1
T(N) = T2
Case 2
T(N) = (2 * a2 * tau * lamda2 * beta(N - 1) - 2 * a2 * tau * h * q2 + h ^ 2 * lamda2 * T(N)) / (lamda2 * h ^ 2 + 2 * a2 * tau * lamda2 * (1 - alfa(N - 1)))
Case 3
T(N) = (lamda2 * h ^ 2 * T(N) + 2 * a2 * tau * (lamda2 * beta(N - 1) + h * k2 * Te2)) / (lamda2 * h ^ 2 + 2 * a2 * tau * (h * k2 + lamda2 * (1 - alfa(N - 1))))
End Select
'Определяем поле температур
For i = (N - 1) To 1 Step -1
T(i) = alfa(i) * T(i + 1) + beta(i)
Next i
Tt(time_i) = T(Ng - 1)
Wend
'Записываем результаты T(i) в таблицу
Dim S As String 'пользовательская переменная
Dim v_Excel As Excel.Application 'область памяти приложения Excel
Dim v_Wb1 As Excel.Workbook 'область рабочей книги
Dim sh As Excel.Worksheet 'область рабочего листа
Set v_Excel = New Excel.Application 'Открываем новое приложение Excel
'Загружаем рабочую книгу
Set v_Wb1 = v_Excel.Workbooks.Open("I:temp.xlsx")
Set sh = v_Excel.Sheets(1) 'Загружаем рабочий лист 1
v_Wb1.Activate 'Получаем доступ к рабочей книге
For i = 1 To N
sh.Cells(i, 1).Value = (i - 1) * h
sh.Cells(i, 2).Value = T(i)
Next i
For i = 1 To 100
sh.Cells(i, 3).Value = (i - 1) * tau
sh.Cells(i, 4).Value = Tt(i)
Next i
v_Excel.Quit 'Закрытие приложения Excel
Set v_Excel = Nothing
End Sub
Список используемых источников
-
Самарский А. А., Вабищевич П. Н. Вычислительная теплопередача. - М.: Едиториал УРСС, 2003. – 782 с. -
Самарский А. А. Теория разностных схем. – М.: Наука, 1977. – 656 с. -
Волков Е. А. Численные методы. – М.: Наука, 1987. – 248 с. -
Вержбицкий В. М. Основы численных методов. – М.: Высшая школа, -
Лыков А. В. Теория теплопроводности. – М.: Высшая школа, 1967. – 600 с.
