Кафедра информационные системы и математическое моделирование

МИНИСТЕРСТВО ОБРАЗОВАНИЯ И НАУКИ РЕСПУБЛИКИ

КАЗАХСТАН

 

МЕЖДУНАРОДНЫЙ УНИВЕРСИТЕТ ИНФОРМАЦИОННЫХ ТЕХНОЛОГИЙ

 

 

КАФЕДРА ИНФОРМАЦИОННЫЕ СИСТЕМЫ И МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ

 

Рысбайулы Б.

 

 

О П О Р Н Ы Й К О Н С П Е К Т Л Е К Ц И Й

По дисциплине Численные методы

.

(Для студентов специальности "математическое и компьютерное моделирование")

 

 

 

АЛМАТЫ 2014

1. Приближенное вычисление определенного интеграла

 

В настоящем пункте рассматриваются способы приближенного вычисления определенных интегралов

 

Введем на [а, в] равномерную сетку с шагом h, т.е. множества точек

 

 

и представим интеграл в виде суммы интегралов по частным отрезкам:

 

 

Для построения формулы численного интегрирования на всем отрезке [а, в] достаточно построить квадратичную формулу для на частном отрезке [х_i-1, х_i].

 

 

Формула прямоугольников

Заменим интеграл S_i выражением Геометрический такая замена означает, что площадь криволинейной трапеции АВСД заменяется площадью прямоугольника АВС₁Д₁ (см. рис. 1).

 

Рис. 1

 

Тогда получим формулу

 

(26)

 

которая называется формулой прямоугольников на частичном отрезке [х_i-1, х_i].

Погрешность метода (26) определяется величиной

 

 

которую легко оценить с помощью формулы Тейлора. Действительно, запишем ψ_i в виде

 

 

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

 

 

Обозначая оценим ψ_i следующим образом:

 

 

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

 

 

т.е. формула имеет погрешность О(h³) при h→0.

Суммируя равенства (26) по I от 1 до N, получим составную формулу прямоугольников

 

 

Погрешность этой формулы

 

 

Отсюда, обозначая получим

 

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

Определение. Приближенное равенство

 

.

 

Называется квадратурной формулой.

 

Формула трапеции

 

На частичном отрезке (х_i-1, x_i) площадь криволинейной трапеции АВСД заменяется площадью прямоугольной трапеции АВСД (рис. 2).

 

 

Рис. 2

 

Тогда

 

Для оценки погрешности

 

Представим его в виде

 

 

Отсюда получим

Составная формула трапеции имеет вид

 

 

Погрешность этой формулы оценивается следующим образом:

 

 

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

Применение формулы трапеции или прямоугольников требует оценки второй производной на отрезке [а, в]. Если такая оценка затруднительна (или вообще невозможно, например, в случае функции определяемых опытным путем), то в предположений малого изменения (или монотонности) второй производной можно во всех полученных оценках заменить множителя М₂h² наибольшей величиной

 

 

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

Для уточнения величины интеграла можно использовать, то обстоятельство, что с уменьшением шага h в два раза погрешность формулы трапеций уменьшается примерно в четыре раза. Отсюда следует, что совпадающие знаки в значениях интеграла, вычисленных с шагом h и можно считать верным. Действительно, если погрешность значения интеграла, вычисленного с шагом обозначить через ε, то погрешность значения интеграла, вычисленного с шагом h, будет приближенно равна 4ε, и значить, разность указанных значений интеграла будет не менее чем 3ε. Поэтому из совпадения m десятичных знаков у рассматриваемых значений интеграла можно заключить, что погрешность , а это означает, что в значений интеграла вычисленном с шагом , все m десятичных знаков верны (здесь предполагается, что погрешность исходных данных пренебрежимо мало).

 

Формула Симпсона

 

При аппроксимации интеграла заменяем функцию f(x) параболой, проходящей через точки (x_I, f(x_I)), I = i-1, i-0,5, i, т.е. представим приближенно f(x) в виде

 

 

Тогда

 

(27)

 

Вычислим

 

 

Из (27) получим, что

 

Таким образом, приходим к приближенному равенству

 

 

которое называется формулой Симпсона.

Погрешность этой формулы ψ_i оценивается так [1]:

 

 

На всем отрезке [a, в] формула Симпсона имеет вид

 

 

Погрешность этой формулы оценивается неравенством:

 

 

Из этой оценки видно, что с уменьшением шага h в два раза погрешность формулы Симпсона уменьшается примерно в 16 раз; поэтому значение интеграла, вычисленное с шагом содержащий на один верный знак больше, чем значение интеграла, вычисленное с шагом h. Это правило на практике очень удобно при оценке точности интеграла.

1.4. Задача 1

Между двумя параллельными сбросами и находится нефтяная залежь В (рис.42) за пределами которой расположены бесконечно простирающая водоносная область. Стрелками показан приток воды из законтурной области. Ширина залежи в = 1000м, толщина пласта h =15м, проницаемость водоносной области k = 0,2·10^-12м², вязкость законтурной воды Упругоемкости β как нефтяной, так и водоносной частей одинаковы, причем β = 2,5·10^-10 Па^-1, вязкость нефти μ_н = 2мПа·С.

 

Рис. 3

 

Отбор жидкости из залежи изменяется во времени следующим образом

 

 

где – время ввода месторождения в разработку. Требуется определить изменение давления на контуре нефтеносности , т.е. по сравнению с начальным давлением после начала разработки залежи.

 

 

Решение. В начале определим пьезопроводность пласта по формуле

 

 

Для расчета изменения во времени давления на контуре нефтяной залежи используя аппроксимацию Карслоу и Егеря [2] имеем:

 

 

 

Данный интеграл вычисления одним из методов: метод прямоугольников, трапеции или Симпсона.

Решение задачи 2

 

 

Рис. 5

 

Для расчета изменения во времени давления на контуре нефтяной залежи, используя аппроксимацию соответствующих решений Карслоу и Егера, имеем

 

 

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

Однако, по условиям данной задачи в период разбурения месторождения объем воды, поступающей из законтурной области, и следовательно, отбираемой жидкости из пласта – переменные во времени. Поэтому для расчета давления на контуре нефтяного месторождения P_кон (t) необходимо использовать интеграл Дюамеля, согласно которому

 

 

В условиях задачи q_ж зависит от физического времени t. В интеграл необходимо поставить Поэтому найдем зависимость q_ж = q_ж(τ) или, что то же самое, q_ж = q_ж(λ). Имеем

 

 

Это формула применяется в том случае, если

 

 

Задание для лабораторной работы.

 

№ п/п	Известные параметры	Определить	Число разбивания
	μ, k, x, H, в, tk, α0, μ, k, x, H, в, tk, α0, μ, k, x, H, в, tk, α0, μ, k, x, H, в, tk, α0, μ, k, x, H, в, tk, α0, μ, k, x, H, R, tk, α0, ρ0, μ, k, x, H, R, tk, α0, ρ0, μ, k, x, H, R, tk, α0, ρ0, μ, k, x, H, R, tk, α0, ρ0,   μ, k, x, H, R, tk, α0, ρ0,	ΔP(t)   ΔP(t)     ΔP(t)   ΔP(t)   ΔP(t)   P(t)   P(t)     P(t)   P(t)     P(t)	n=100, 200   n=200, 400     n=300, 600   n=150, 300   n=250, 500   n=300, 600   n=250, 500     n=250, 500   n=200, 400     n=500, 1000

Структурная схема расчета.

 

 

- - - - - - - - - - ввод начальных данных

t = 0, t_max, Δt ΔP конец

 

Рис.6

Фильтрация жидкости и газа

Теория фильтрации изучает движение газов, жидкостей и их смесей в порах, т.е. в твердых телах, пронизанных системой сообщающих между собой пустот (пор), что делает их проницаемой для жидкостей. Пористая среда состоит огромного числа случайно расположенных зерен различной формы и величины. Поэтому пространство, в котором движется жидкость, представляет собой систему пор, непрерывно переходящих одна в другую. Для пористой среды характерно свойство сообщаемости пор. Из-за нерегулярности строения порового пространства его нельзя полностью описать никаким конечным набором параметров; для целей теории фильтрации, однако, достаточно небольшого числа осредненных характеристик. Важнейшая характеристика пористой среды – ее пористость , равная отношению объема, занятого в выделенном элементе порами, к общему объему элемента .

Самым простым и наиболее изученным случаем нестационарной фильтрации является фильтрация слабосжимаемой жидкости в упругодеформируемом пласте (в технических приложениях эти задачи получили название задач упругого режима фильтрации). В основу исследования кладется система уравнений закона фильтраций и уравнения неразрывности:

 

 

Здесь -вектор скорости фильтрации; - плотность жидкости. Для того, чтобы получить замкнутую систему уравнений, нужно воспользоваться тем, что свойства жидкости (плотность ρ и вязкость μ), так же как и пористость и проницаемость пористой среды, являются функциями давления (мы предположим движение изотермическим). Исходя, из предположения, слабой сжимаемости жидкости и пористой среды получается уравнение [14-15]:

 

где коэффициент .

 

Здесь - коэффициент пьезопроводности. Полученное уравнение обычно называется уравнением упругого режима или уравнением пьезопроводности.

Рассмотрим постановку основных задач теории упругого режима. Определим распределение давления в некоторой замкнутой области пространства D на протяжении промежутка времени 0 ≤ t ≤ T. Из теории уравнения теплопроводности известно, что если задать на границе Г области D линейную комбинацию давления и его производной по нормали к границе области

 

и задать начальное распределение давления в области Д

 

,

 

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

Рассмотрим подробнее физический смысл тех или иных дополнительных условий. Область, в которой ищется распределения давления жидкости, обычно представляет собой пористый пласт, частично имеющий непроницаемые границы, а частично сообщающийся с другими пластами и вскрывающими его скважинами. На непроницаемых границах должно удовлетворяться очевидное условие отсутствие потока – равенство нормального компонента скорости фильтрации нулю:

 

U_n = 0,

 

Откуда используя закон Дарси, получаем.

 

 

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

 

.

 

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

При рассмотрении нестационарных процессов в залежи давление в водоносной области можно считать постоянным. Следует, однако, отчетливо представлять себе, что понятие области постоянного давления не является абсолютным. Чем более длительный характер носят изменения давления, тем на большую область они распространяются.

Часть границы области фильтрации обычно образована стенками скважины и дренажных галерей. На этой части границы чаще всего задается либо давление жидкости, либо поток ее через стенки скважины. Выбор того или иного условия зависит от режима работы скважины или галерей. Могут быть и более сложные условия, когда задается связь давления с расходом жидкости. Задание потока жидкости согласно закону Дарси эквивалентна заданию нормальной производной от давления. Условия этого типа выполняются на тех участках границы, через которые может происходить обмен жидкости с соседними пластами через сравнительно слабопроницаемые перемычки. Если толщина перемычки Δ мала, а давление за ней можно считать постоянным, то расход вытекающей жидкости через участок перемычки площадью ds составит . Это количество жидкости должно быть равно

 

uде U_n – нормальная проекция скорости фильтрации на рассматриваемом участке границы. Отсюда имеем

 

 

т.е. условия третьего рода.

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

Задача 3. Одинарное прямолинейно-параллельное движение. Пусть скорость фильтрации жидкости параллельны оси х и не зависит от координат y и z. Давление при этом удовлетворяет уравнению Предположим, что в плоскости х = L давление сохраняет постоянное значение, равное контурному.

 

При х = 0 задан расход жидкости q(t). При заданных значениях μ, k, x и L определить давление пласта.

Решение. В работе [5], предполагая, что получено решение поставленной задачи. Она описывается в виде:

 

 

Полученная формула имеет двоякий смысл. С одной стороны, он описывает распределение давления в пласте конечной длины L при малом времени. С другой стороны он дает распределение давления в произвольный момент времени. В пласте «бесконечной» протяженностью L → ∞. Дело в том, что конечное (не бесконечно малое) изменения давления распространяется за заданное время лишь на конечное расстояние и, если рассматриваются малые времена, можно считать пласт бесконечным. Решение задачи для бесконечного пласта автомодельное: независимые переменные х и t входят в решение не порознь, а лишь в комбинации

Если , то [5]

 

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

 

P ₁ (0, t) = P₀, P(L, t) = P(L, 0) = const;

 

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

Оценим это время для систем различных режимов. При этом для χ характерное значение х = 10⁴см²/сек. В результате имеем: при L = 1м (переходный процесс в одном пуске-блике породы) t = 0,1сек; при L = 300м (порядка расстояния между скважинами) t = 10⁴сек»3г. L = 1км (порядка размеров месторождения) t = 10⁷сек £ 100 суткам; L = 100км (порядка размеров крупной водонапорной системы) t = 10⁹сек = 30 лет.

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

Задания для лабораторной работы.

 

№ п/п	Известные параметры	Определить	Число разбивания
	μ, k, χ, q(t) = α0 t, t £ T, P0 (x) μ, k, χ, q(t) = α0tsin , t £ T, P0 (x) μ, k, χ, q(t) = α0ttg ,t £ T, P0 (x) μ, k, χ, q(t) = α0t^2/(t^0 + 2) t £ T, P0 (x) μ, х, q(t) = α0t, t £ T, P (t1), P1 (t1+ Δt) k, х, q(t) = α0tsin P (t1), P (t1+ Δt) k, х, q(t) = α0t, ρ (t1), P (t1+ Δt) μ, х, q(t) = α0tsin , ρP(t1), P (t1+ Δt) μ, k, q(t) = α0t, ρ (t1), P(t1+ Δt) μ, k, q(t) = α0tsin , P(t1), P(t1+ Δt)	P (x, t) P (x, t)   P (x, t)   P (x, t) k, ΔP (x, t)   μ, ΔP (x, t)   μ, ΔP (x, t) k, ΔP (x, t)   x, ΔP (x, t) x, ΔP (x, t)	n=200, 400 n=300, 600   n=400, 800   n=500, 1000 n=250, 500   n=350, 700   n=400, 800 n=500, 1000   n=500, 1000 n=600, 1200

Задача 4. Рассмотрим теперь одномерное осесимметричное (плоско-радиальное) движение при упругом режиме. Распределение давления определяется при этом как решение уравнения теплопроводности в полярных координатах (r, j):

 

 

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

 

P (r, 0) = f(r)

 

а граничным условиям при r = ρ и r = R.

По прежнему основной интерес представляют решения, отвечающие стационарному начальному распределению давления

 

f(r)=C₁lnr+C₂.

 

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

 

 

 

где – постоянная размерность давления.

 

 

Рис. 7

 

Предполагая, что в работе [5] получена решение поставленной задачи в виде

(28)

 

Отметим важное обстоятельство: в соотношения для P (r, t) не входит радиус скважины . Это означает, что в области применяемости условия распределения не зависит от радиуса скважины.

(28) – является несобственным интегралом, с помощью замены переменной приведем его собственному интегралу. Для этого сделаем замену:

 

 

Первоначальное дифференциальное уравнение является линейным, поэтому здесь применим принцип суперпозиции. Если имеется две скважины расположенные друг от друга на расстоянии l, то давление пласта в точке А удаленной от первой скважины на расстоянии будет равен (рис. 7)

 

 

Задания для лабораторной работы

№ п/п	Известные параметры	Определить	Число разбивания
	k, χ, μн, q1, q2 , R1, R2, t £ T k, χ, μн, q1, q2 , 0 £ х^2+ у^2 £ R^2, t £ T k, χ, μн, q1, q2 , q3, 0 £ х^2+ у^2 £ R^2, t £ T k, χ, μн, q1, q2 , q3, q4, 0 £ х^2+ у^2 £ R^2, t £ T k, χ, μн, q1=q2, 0 £ х^2+ у^2 £ R^2, t £ T k, χ, μн, q1= q2 = q3, 0 £ х^2+ у^2 £ R^2, t £ T	ΔP (t) ΔP(x, у, t) ΔP(x, у, t) ΔP(x, у, t) ΔP(x, у, t) ΔP(x, у, t)	n=500 n=300 n=200 n=300 n=250 n=400

 

Сплайн 2-го порядка S(x)

На каждом из отрезков (x_i, x_{i +1}) функция у = f(x) приближается параболой

 

S(x) =

В узлах х = х_i ставятся следующие условия

 

1)

 

2) - непрерывность функции S (x) в узлах

 

x_i , i = 1, 2, …, n-1.

3) - непрерывность первой производной функции

S (x) в узлах x_i , i = 1, 2, …, n-1.

 

Используя 1) - 3) определяем а_i, b_i, с_i.

 

1. Из условия следует, что

.

Сюда можно дописать дополнительный коэффициент .

2. Из условия 2) и 3) получаем систему

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

Задача. Длинный трубопровод с теплопроводностью λ (ккал/м.час град.) находится в состоянии теплового равновесия, т.е. температура точек трубопровода не изменяется во времени. Потеря тепла через поверхность трубопровода в окружающую среду, температура которой , пропорциональна разности температур с постоянным коэффициентом теплопередачи α (ккал/м². час град.). Считая температуру θ во всех точах поперечного сечения трубопровода постоянной, найти ее зависимость θ = θ(х) от координаты, отсчитываемой от какого-либо конца.

 

Математическая модель

 

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

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

Пусть длина стержня равна l м, периметр поперечного сечения p м, площадь поперечного сечения Q м². Исследуем процесс распространения тепла в элементарном отрезке длиной dx на расстоянии х от того конца стержня, температура которого t₁. Количество тепла проходящего за время dt через сечение трубопровода находящееся на расстоянии х от начало трубопровода, согласно теории теплопередачи, будет равно:

 

 

Количество тепла, прошедшее за время dt через сечение, находящееся на расстоянии х + dx от начала, будет равна:

 

 

Участок стержня, заключенный между сечениями, отстоящими от начала на расстояниях х и х + dx, вследствие теплопроводности, приобретает за время dt количество тепла, равное разности указанных количеств, т.е.:

 

 

За то же время потеря тепла от этого же участка в окружающую среду будет равна:

 

 

Но так как изучаемый нами процесс является стационарным, то:

 

 

Окончательно мы приходим к следующему дифференциальному уравнению:

 

 

В итоге получена задача:

 

Это уравнение с постоянными коэффициентами. Его общий интеграл будет:

 

 

Используя граничные условия, составим систему:

 

откуда

 

Подставляя значения С₁ и С₂ получим:

 

 

Выделим элемент длины dх, находящийся на расстояний х от левого конца, и примем его температуру равной . За время ∆t через левую границу этого элемента пройдет количество тепла

а через правую на расстоянии х+dх от конца

Таким образом, выделенный участок приобретает за время t количество тепла, равное разности

.

Вместе с тем потеря тепла этого элемента через поверхности в окружающую среду равна

.

Вследствие предположения о стационарности процесса эти количества равны, т.е.

откуда

(4.1)

Пусть на обоих концах стержня поддерживается постоянная температура . Тогда краевые условия имеют вид.

, (4.2)

 

Численный пример. Пусть a = 10

l = 300 ккал/м×час×град

 

тогда

 

При этом случае получится зависимость

 

 

Переменные. Блок-схема

 

Из определения А_і, В_і, С_і следует равенства

 

С_і = А_{і -1}, В_і = А_і +А_{і - 1} + а² h²

Поэтому, (4.9) можно переписать в виде

 

На основе этих формул, при программировании используются массивы

A[0…n], α[0…n-1], β[0…n-1].

БЛОК-СХЕМА

 

αn-l =0, βn-1 = θ2

 

l=n-1, 0,-1

		конец

 

 

l=0, n- 1, 1

 

Цель лабораторной работы.

 

С помощью программы установить:

- если длина трубы l достаточно большой, то граничные условия практически не влияют на распределение температуры вдоль трубопровода.

- если длина трубы l короткая, то θ₁ и θ₂ влияет на распределение температуры вдоль трубопровода.

- изучить влияние α и θ₀ на распределение температуры.

 

 

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

Подогретый до температуры θ₁ (град.) нефть перевозится по подземному трубопроводу.

Глубина прокладки нефтепровода Н(м).

 

Считая, что теплоотдача происходит только по вертикальному направлению определить распределение температуры от трубы до поверхности земли.

х     Н   О

Рис. 12

Математическая модель.

Ось ОХ направ