Трушков А.С. Учебно-методический комплекс по дисциплине Компьютерное моделирование

Подождите немного. Документ загружается.

б) Для модели из пункта а) исследовать зависимость температуры на ле-

вой границе стержня от интенсивности теплоисточника f и его длины

c = p

- l

, принимая в качестве закона теплообмена на левой границе

′

′′

′

(0) = 0,9(Т - 30

). Построить таблицы и графики зависимости u(0, f, c), привес-

ти характерные диаграммы развития процесса.

Пример решения задачи.

Задание. Исследовать задачу распространения тепла в стержне с тепло-

изолированной боковой поверхностью.

Будем решать задачу в табличном процессоре MS Excel. Оформим рабо-

чий лист Расчет для ввода исходных данных и оформим строку заголовков для

таблицы результатов, как показано на рис. 1. Вычисления организуем в макросе

Visual Basic for Applications (VBA), текст которого с комментариями приведен

ниже.

' Инструкция об обязательном объявлении переменных

Option Explicit

' Объявление объектовых переменных

' ID – массив исходных данных

' Rez – массив результатов

' Graf – массив исходных данных для динамической графики

Dim ID As Object

Dim Rez As Object

Dim Graf As Object

' Объявление переменных

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

' lmax, Tmax – длина стержня и время счета

' dl, dt – шаги разностной сетки по длине и времени

' a0, a1, b0, b1 – коэффициенты левой части краевых условий

' nlv, ntv - количество шагов промежуточного вывода по длине и времени

' AA, BB - коэффициенты правой части краевых условий

' lg, pg, fip – левая и правая границы теплоисточника, его интенсивность

' x, t – текущие значения длины и времени

' nl, nt – число шагов разностной сетки по длине и времени

' priz – признак расчета, priz = 1 – динамическая графика

Dim a As Double, lmax As Double, Tmax As Double

Dim dl As Double, dt As Double

Dim a0 As Double, a1 As Double, b0 As Double, b1 As Dou-

ble

Dim nlv As Integer, ntv As Integer, AA As Double

Dim BB As Double

Dim lg As Double, pg As Double, fip As Double, x As Dou-

ble

Dim t As Double, nl As Integer, nt As Integer

Dim priz As Integer

' Процедура обработки события Click кнопки Пуск

Sub Пуск()

priz = 0

Расчет

End Sub

' Процедура обработки события Click кнопки Графика

Sub Графика()

priz = 1

Расчет

End Sub

' Функция задания интенсивности теплоисточника f

Function fii(i)

x = i * dl

If (x > lg) And (x < pg) Then fii = fip Else fii = 0

End Function

' Функция задания начального распределения температуры

Function u0(i)

u0 = 0

End Function

' Расчетная процедура

Sub Расчет()

' Объявление переменных

' k1, k2, m1, m2 – коэффициенты разностной схемы (10) для краевых условий

' j, i – переменные цикла

' ktv – текущая строка вывода результатов

' t_name, t_range – имя и диапазон исходных данных диаграммы динамической

' графики

Dim k1 As Double, k2 As Double, m1 As Double, m2 As Dou-

ble

Dim j As Integer, i As Integer, ktv As Integer

Dim t_name, t_range

' Объявление динамических массивов для коэффициентов разностной схемы

(8)-(9)

Dim ai() As Double

Dim bi() As Double

Dim ci() As Double

Dim Fi() As Double

' Объявление динамического массива для искомой температуры

Dim y() As Double

' Объявление динамических массивов прогоночных коэффициентов (13)-(15)

Dim alf() As Double

Dim bet() As Double

' Активация листа Расчет

Sheets("Расчет").Activate

' Связывание объектовых переменных с диапазонами данных

Set ID = Worksheets("Расчет").Range("A3:A15")

Set Rez = Worksheets("Расчет").Range("A24:L20000")

Set Graf = Worksheets("Графика").Range("A1:B20000")

' Обнуление содержания объектовой переменной Rez

Rez.ClearContents

' Присвоение значений исходным данным из объектовой переменной ID

a = ID(1)

lmax = ID(2)

Tmax = ID(3)

nl = ID(4)

nt = ID(5)

a0 = ID(6)

a1 = ID(7)

b0 = ID(8)

b1 = ID(9)

AA = ID(10)

BB = ID(11)

nlv = ID(12)

ntv = ID(13)

lg = ID(14)

pg = ID(15)

fip = ID(16)

' Вычисление шагов разностной сетки

dl = lmax / nl

dt = Tmax / nt

' Определение размерности динамических массивов

ReDim ai(1 To nl - 1)

ReDim bi(1 To nl - 1)

ReDim ci(1 To nl - 1)

ReDim Fi(1 To nl - 1)

ReDim y(nl)

ReDim alf(1 To nl)

ReDim bet(1 To nl)

' Вывод результатов в начальный момент времени с учетом частичного выво-

да

' по длине стержня

Rez(1, 1) = 0

Rez(1, 2) = u0(0)

For i = 1 To nl / nlv

Rez(1, i + 2) = u0(i * nlv)

Next i

' Графическая обработка в начальный момент времени

If priz = 1 Then

' Активация листа Графика

Sheets("Графика").Activate

' Обнуление содержания объектовой переменной Graf

Graf.ClearContents

Graf(1, 2) = "t=0"

' Формирование исходных данных для динамической графики в начальный мо-

мент

For i = 0 To nl

Graf(i + 2, 1) = i * dl

Graf(i + 2, 2) = u0(i)

Next i

' Отработка динамической графики в начальный момент

t_name = ActiveSheet.ChartObjects(1).Name

t_range = "A1:B" & CStr(1 + nl)

ActiveSheet.ChartObjects(t_name).Activate

ActiveChart.ChartArea.Select

ActiveChart.SetSourceData

Source:=Sheets("Графика").Range(t_range), PlotBy _

:=xlColumns

End If

' Определение начальных значений

ktv = 1

t = 0

j = 1

For i = 0 To nl

y(i) = u0(i)

Next i

' Цикл реализации метода прогонки

Do Until t > Tmax

' Текщее время

t = t + dt

' Вычисление коэффициентов разностной схемы (8)-(9)

For i = 1 To nl - 1

ai(i) = a ^ 2 * dt / 2 / dl ^ 2

bi(i) = a ^ 2 * dt / 2 / dl ^ 2

ci(i) = 1 + a ^ 2 * dt / dl ^ 2

Fi(i) = (1 - a ^ 2 * dt / dl ^ 2) * y(i) + _

a ^ 2 * dt / 2 / dl ^ 2 * (y(i - 1) + y(i + 1))

+ dt * fii(i)

Next i

' Вычисление коэффициентов (10) разностной схемы для краевых условий

k1 = -a1 / (a0 * dl - a1)

m1 = AA * dl / (a0 * dl - a1)

k2 = b1 / (b0 * dl + b1)

m2 = BB * dl / (b0 * dl + b1)

' Прямая прогонка (12)-(15)

alf(1) = k1

bet(1) = m1

For i = 1 To nl - 1

alf(i + 1) = bi(i) / (ci(i) - ai(i) * alf(i))

bet(i + 1) = (ai(i) * bet(i) + Fi(i)) / (ci(i) -

ai(i) * alf(i))

Next i

' Обратная прогонка (16)-(17)

y(nl) = (m2 + k2 * bet(nl)) / (1 - alf(nl) * k2)

For i = nl - 1 To 0 Step -1

y(i) = alf(i + 1) * y(i + 1) + bet(i + 1)

Next i

' Вывод результатов в текущий момент времени с учетом частичного вывода

' по длине стержня и времени

If j = ntv Then

ktv = ktv + 1

Rez(ktv, 1) = t

Rez(ktv, 2) = y(0)

For i = 1 To nl / nlv

Rez(ktv, i + 2) = y(i * nlv)

Next i

j = 1

Else

j = j + 1

End If

' Графическая обработка в текущий момент времени

If priz = 1 Then

' Формирование исходных данных для динамической графики в текущий мо-

мент

Graf(1, 2) = "t=" + CStr(Format(t, "0.00"))

For i = 0 To nl

Graf(i + 2, 2) = y(i)

Next i

' Отработка динамической графики в текущий момент

ActiveSheet.ChartObjects(t_name).Activate

ActiveChart.ChartArea.Select

ActiveChart.SetSourceData

Source:=Sheets("Графика").Range(t_range), PlotBy _

:=xlColumns

End If

Loop

End Sub

Обработка результатов после проведения расчета осуществляется при

нажатии кнопки Пуск, созданной с помощью панели инструментов Формы и

связанной с процедурой Пуск. До начала расчета создается и форматируется

диаграмма, исходные данные которой соответствуют выбранным моментам

времени из массива результатов на рабочем листе Расчет (рис. 2).

В программе предусмотрено динамическое объявление массивов, размер-

ности которых определяются после ввода размерностей разностной сетки.

Реализация решения разностной системы уравнений с использованием

разностной схемы Кранка-Николсона осуществляется в процедуре Расчет.

Исследование точности разностной схемы. Рассмотрим исследуемую

модель при следующих исходных данных:

a = 0,5; l = 1; T = 1; f(x,t) = 0; u

(x) = 0;

= 1; a

= 0; A = 30; b

= 1; b

= 0; B = 0.

Оценим по правилу Рунге погрешность разностной схемы Кранка-

Николсона в каждой точке x

стержня, учитывая второй порядок аппроксима-

ции по пространственной производной:

)h(y)2/h(y

)2/h(y)x(yd

−

−−

−

−−

−

−−

−

≈

≈≈

≈−

−−

−=

Будем считать, что на разностной сетке большой размерности прибли-

женное решение практически совпадает с точным. В качестве такой сетки вы-

берем N = 2000 – число шагов по пространственной переменной, L = 2000 –

число шагов по времени. Для получения приближенного решения y

(h)

рас-

смотрим сетку N

××

L = 100

××

2000, а приближенного решения y

(h/2) -

××

L =

200

××

2000. В качестве оценки погрешности выбирается максимальная из точеч-

ных погрешностей во всем исследуемом интервале времени. В качестве факти-

ческой погрешности принимается максимальное рассогласование результатов,

полученных на сетке 200

××

2000 и 2000

××

2000. Результаты расчетов на указан-

ных сетках, а также вычисления оценок погрешностей по правилу Рунге и их

фактических величин в выбранных узлах разностной сетки размещаются в ра-

бочем листе Правило Рунге и приведены на рис 3. В ячейке N34 находится

оценка максимальной погрешности по правилу Рунге, а в ячейке N64 – макси-

мальная фактическая погрешность. Для рассматриваемого примера оценка по-

грешности составляет d = 0,000398, а ее фактическое значение d

= 0,000394.

Хорошее соответствие этих результатов подтверждает второй порядок аппрок-

симации разностной схемы Кранка-Николсона.

Исследование процесса теплопроводности в стержне. Рассмотрим раз-

личные процессы теплопроводности в стержне.

а) Рассмотрим модель при следующих исходных данных:

a = 0,5; l = 1; T = 1; f(x,t) = 0; u

(x) = 0;

= 1; a

= 0; A = 40; b

= 1; b

= 0; B = 0.

Это значит, что на левой границе задана температура u(0) = 40, на правой

границе u(1) = 0, теплоисточники в стержне отсутствуют, начальная температу-

ра u

(x) = 0. Графическая обработка процесса теплообмена приведена на рис. 4.

б) Рассмотрим модель при следующих исходных данных:

a = 0,5; l = 1; T = 1; f(x,t) = 0; u

(x) = 0;

= 1; a

= 0; A = 40; b

= 1; b

= 1; B = 10.

По сравнению с моделью а) на правой границе задано условие теплооб-

мена в виде:

′

′′

′

(1) = 10 - T.

Графическая обработка процесса теплообмена приведена на рис. 5.

с) Рассмотрим модель при следующих исходных данных:

a = 0,5; l = 1; T = 1; u

(x) = 0;

[

[[

[

]

]]

]

[

[[

[ ]

]]

]



















∉

∉∉

∉

∈

∈∈

∈

.7,0;5,0x,0

;7,0;5,0x,100

)t,x(f

= 1; a

= 0; A = 40; b

= 1; b

= 1; B = 10.

По сравнению с моделью б) внутри стержня задан теплоисточник с ин-

тенсивностью f = 100 (т.е. внутренний теплоподвод) на отрезке с левой грани-

цей l

= 0,5 и правой границей p

= 0,7. Графическая обработка процесса тепло-

обмена приведена на рис. 6.

Рис. 1.

Рис. 2.

Рис. 3.

Рис. 4.

Рис. 5.