Температурный расчет с помощью вычислений информационной математики

Постановка задачи.

По длинной квадратного сечения трубе течет горячая жидкость. Труба наполовину погружена в ледяную ванну, так, что температура нижней половины поверхности трубы равна 0⁰ С. Верхняя плоскость трубы имеет постоянную температуру 100 ⁰ С. На участке между ледяной ванной и верхней плоскостью температура наружной поверхности трубы изменяется линейно по высоте от 0 ⁰ С до 100 ⁰ С. Жидкость внутри трубы имеет температуру 200 ⁰ С.

Рис. 3.

Распределение температуры в теле трубы удовлетворяет уравнению

С погрешностью не более 0,5 ⁰ С вычислить распределение температуры в теле трубы.

Математическая формулировка задачи.

Решить диф.уравнение в частных производных:

с задаными началиными условиями на границах области дифференцирования.

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

в результате чего диф.уравнение преобразуется в 5-ти диаганальную систему алгеброических уравнений n-го порядка.

Систему алгеброических уравнений буду решать методом Зейделя.

Погрешность решения задачи найду по формуле:

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

Функциональная схема.

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

Описание метода.

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

Конечные разности и производные.Пусть некоторая функция y(x) задана на отрезке [a,b]. Будем считать, что она непрерывна и многократно дифференциру­ема на этом отрезке. Разделим отрезок на равные части длиною h и обозначим точки деления x0,x1,...,xi,...,xn.Значе­ния функции в этих точках обозначим соответственно y0,y1,...,yi,...,yn.Первой центральной разностью в i-й точке (i=1,2,...,n-1) называют разность:

С помощью этой разности можно приближенно вычислить значение первой производной у` в i-й точке.

Разложим функцию y(x) в степенной ряд. приняв за центр разложения точку xi и ограничившись четырьмя членами:

где

Аналогично найдем значение ф-ции и в точке,отстоящей от центра разложения на шаг (-h):

где .

Подставляя получим:

Таким образом,производная y` приближонно заменяется конечно-разностным отношением с ошибкой порядка h*h:

Второй центральной разностью ф-ции y(x) в i-й точке называют величину:

С помощью этой разности можно приближонно вычислить значение второй производной y`` в i-й точке.Используем теперь 5 членов разложения в ряд Тейлора:

Таким образом,вторая производная y`` с ошибкой порядка h*h может быть приближонно заменена конечно-разностным отношением:

При определении разностей в i -и точке использовались значения функции в точках, расположенных симметрично относительно xi . Поэтому эти разности назы­ваются центральными.

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

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

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

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

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

4.В разностной форме записывают краевые условия и состав­ляют полную систему разностных уравнений.

Оценка погрешности решения краевой задачи

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

Для линейных краевых задач доказана теорема о том, что по­рядок точности решения краевой задачи не ниже порядка точности аппроксимации производных конечно-разностными отношениями. Оценку погрешности производят при­емом Рунге. Краевую задачу решают дважды: с шагом сетки h и с шагом сетки H=kh, погрешность решения с малым шагом h оценивают по формуле:

где y(h) и y(H) - решения, полученные для одной и той же точ­ки -xi отрезка интегрирования с разными шагами. Относительную погрешность E оценивают по формуле:

Если при составлении разностной системы уравнений исполь­зуются левые или правые разности, то погрешность решения будет выше, порядка 0(h), и для ее оценки в формулах следует заменить k*k на k .

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

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

Блок схема.

Подпрограмма МКР.

c------------------------------------------------------------------

c ПОДПРОГРАММА СОСТАВЛЕНИЯ СИСТЕМЫ УРАВНЕНИЙ

c МЕТОДОМ КОНЕЧНЫХ РАЗНОСТЕЙ

c

c real H-шаг по оси X

c real K-шаг по оси Y

c real N-количество уравнений(примерное число,желательно N=M*P)

c real y(6,N)-выходной массив уравнений,содержащий следующие поля:

c y(1,N)-номер точки по оси X

c y(2,N)-номер точки по оси Y

c y(3,N)-коэфициен уравнения для Q(y(1,N)-1,y(2,N))

c y(3,N)=h^2/(2*(h^2+k^2))

c y(4,N)-коэфициен уравнения для Q(y(1,N),y(2,N)-1)

c y(4,N)=k^2/(2*(h^2+k^2))

c y(5,N)-коэфициен уравнения для Q(y(1,N)+1,y(2,N))

c y(5,N)=h^2/(2*(h^2+k^2))

c y(6,N)-коэфициен уравнения для Q(y(1,N),y(2,N)+1)

c y(6,N)=k^2/(2*(h^2+k^2))

c integer M-число узлов по оси X

c integer P-число узлов по оси Y

c real Q(M,P)-массив значений Y

c integer N-выходное количество получившихся уравнений

c------------------------------------------------------------------

subroutine mkr(H,K,N,y,M,P,q)

integer M,P,IIX,IIY,NN,N,KR1,KR2,KR3

real y(6,N),H,K,q(M,P),HX,KY

c-----------------------------------------------------------------

c подсчитываю коэфициенты

c h^2/(2*(h^2+k^2))

c и

c k^2/(2*(h^2+k^2))

c-----------------------------------------------------------------

HX=H**2/(2*(H**2+K**2))

KY=K**2/(2*(H**2+K**2))

c-----------------------------------------------------------------

c составление уравнений

c и

c присваивание начальных значений

c

c nn-счетчик уровнений

c iix-номер текущего узла по оси X

c iiy-номер текущего узла по оси Y

c-----------------------------------------------------------------

NN=0

KR1=((P-1)/8)*3+1

KR2=((P-1)/8)*5+1

KR3=((M-1)/4)*3+1

do IIY=2,P-1

do IIX=2,M

if (NN.eq.N)then

print *,'ПЕРЕПОЛНЕНИЕ МАССИВА Y'

stop

endif

c-----------------------------------------------------------------

c проверка границы трубы с жидкостью

c-----------------------------------------------------------------

if ((IIY.ge.KR1).and.(IIY.le.KR2).and.(IIX.ge.KR3)) then

q(IIX,IIY)=200.

c-----------------------------------------------------------------

c