Форум программистов, компьютерный форум, киберфорум
Наши страницы
WH
Войти
Регистрация
Восстановить пароль
Оценить эту запись

Истечение жидкости из резервуара

Запись от WH размещена 03.01.2020 в 11:41
Обновил(-а) WH Вчера в 04:18

Очень часто в интернете встречается задача на определение времени истечения жидкости из резервуара. Для примера была взята задача приведенная на этом форуме, истечение жидкости из шарового резервуара с объемом 600 м.куб и отверстием в днище диаметром 0.12 м. По умолчанию будем полагать, что сверху резервуар имеет "дыхательное" отверстие достаточного сечения. Для решения задачи нужно решить дифференциальное уравнение:

http://www.cyberforum.ru/cgi-bin/latex.cgi?\frac{dx}{dt} = - \frac{{S}_{o} \cdot K \cdot \sqrt{2gx}}{S(x)}

где:

http://www.cyberforum.ru/cgi-bin/latex.cgi?S(x) - площадь сечения шара на высоте http://www.cyberforum.ru/cgi-bin/latex.cgi?x,
http://www.cyberforum.ru/cgi-bin/latex.cgi?{S}_{o} - сечение отверстия в днище,
http://www.cyberforum.ru/cgi-bin/latex.cgi?K - коэффициент скорости истечения жидкости через насадки (примем равным 1),
http://www.cyberforum.ru/cgi-bin/latex.cgi?x - высота до поверхности жидкости,
http://www.cyberforum.ru/cgi-bin/latex.cgi?t - время,
http://www.cyberforum.ru/cgi-bin/latex.cgi?g - ускорение свободного падения.

Коэффициент на вязкость жидкости учитывать не будем, поскольку задача учебная. Подкоренное выражение http://www.cyberforum.ru/cgi-bin/latex.cgi?\sqrt{2gx} представляет собой формулу Торичелли, это скорость истечения жидкости. Знак минус означает, что уровень жидкости уменьшается по мере вытекания жидкости из резервуара.

По сути нам нужно решить задачу Коши, решим ее численным способом с использованием самого простого метода решения диф. уравнений - метода Эйлера. Код решения на языке Ada.

код на языке Ada

Pascal
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
--Определение времени слива жидкости из шарового резервуара
with Ada.Text_IO; use Ada.Text_IO;
with Ada.Float_Text_IO; use  Ada.Float_Text_IO;
with Ada.Numerics.Elementary_Functions; use Ada.Numerics.Elementary_Functions;
--with Ada.Integer_Text_IO; use Ada.Integer_Text_IO;
 
procedure Main is
 
--объявим переменные
 pi                     :constant := 3.1416;    --число ПИ
 g                      :constant := 9.8067;    --ускорение свободного падения
 dotv, sotv             :float;                 --Диаметр и сечение отверстия в резервуаре
 vshar, rshar, dshar    :float;                 --Объем и радиус и диаметр шарового резервуара
 z, dz                  :float;                 --Высота поверхности жидкости z и ее приращение
 t, dt                  :float;                 --Время t и его приращение
 s                      :float;                 --Площадь поперечного сечения резервуара на высоте z
 k                      :float;                 --Коэффициент скорости истечения жидкости
 h                      :float;                 --Промежуточное значение для вычисления s 
 
--расчет
begin        
 
  --Начальные параметры резервуара и отверстия
    vshar:= 600.0; dotv:= 0.12;
    
  --Расчет исходных параметров резервуара
    sotv:= pi * dotv * dotv / 4.0;
    rshar:= (3.0 * vshar / (4.0 * pi))**(1.0/3.0);
    dshar:= rshar * 2.0;
    z:= rshar * 2.0;
    
  --Исходные значения для решения уравнения
    t := 0.0; dz := 1.0e-3; k :=1.0;
     
  --Проинтегрируем диф. уравнение   
    while z > dz loop
      h  := z-rshar;
      s  := pi * (rshar**2 - h**2);
      dt := dz * s / (sotv * k * sqrt(2.0 * g * z));
      t  := t + dt; 
      z  := z - dz;
    end loop;
       
 --вывод результата      
 Put ("Объем резервуара         = "); Put (vshar); Put (" м.куб.");  New_Line;
 Put ("Диаметр резервуара       = "); Put (dshar); Put (" м.");      New_Line;
 Put ("Сечение отверстия        = "); Put (sotv ); Put (" м.кв.");   New_Line;
 Put ("Время вытекания жидкости = "); Put (t);Put (" с.");
 
end Main;

Результат:
Код:
Объем резервуара         =  6.00000E+02 м.куб.
Диаметр резервуара       =  1.04645E+01 м.
Сечение отверстия        =  1.13098E-02 м.кв.
Время вытекания жидкости =  5.92495E+03 с.
Шаг интегрирования в данной программе задан - 1 мм по высоте уровня жидкости, это значит, что общее число шагов интегрирования получилось 10464 (диаметр резервуара в мм) / 1 = 10464. Приведенную программу можно запустить в том числе с использованием online компилятора и поэкспериментировать с разными параметрами. Запуск осуществляется нажатием на кнопку Run.

Ниже представлено то же самое решение на языке Фортран. Однако переписывать один к одному не столь интересно, поэтому переписано с применением объектно-ориентированной технологии. Имеем объект - шар с набором присущих ему параметров и к нему метод вычисления времени истечения жидкости. Создали объект, задали ему параметры, вызвали метод вычисления времени и получили ответ. В этом же классе можно прописать параметры и методы для резервуаров любой другой формы. Здесь правильно будет сказать, что и на языке Ada можно писать программы с применением ООП.

модуль с классом

Fortran
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
module telo
implicit none
public :: shar, method_t                        !Делаем класс и его методы публичными
private                                         !При этом сам класс будет черным ящиком
 
    type :: shar
     real       :: dotv                         !Диаметр отверстия в резервуаре
     real       :: vshar                        !Объем резервуара
     real       :: k = 1.0                      !Коэффициент скорости истечения жидкости
     real       :: t = 0.0                      !Время истечения жидкости
    
    contains                                    !Методы класса
     procedure :: method_t => method_t          !вычисление времени истечения жидкости
    end type shar
 
 contains
    
  subroutine method_t (this)
  class (shar) :: this
    real, parameter :: pi=3.1416, g = 9.8067    !Константы: число "Пи" и ускорение свободного падения "g"
    real        :: sotv                         !Cечение отверстия в резервуаре
    real        :: rshar                        !Радиус шарового резервуара
    real        :: z, dz                        !Высота поверхности жидкости z и ее приращение
    real        :: dt                           !приращение времени
    real        :: s                            !Площадь поперечного сечения резервуара на высоте z
    
    !Расчет   исходных параметров резервуара
    sotv    = pi * this%dotv * this%dotv / 4.0
    rshar   = (3 * this%vshar / (4 * pi)) ** (1.0/3.0)
    z   = rshar * 2.0
 
    !Исходные значения для решения уравнения
    this%t = 0.0; dz = 1e-3
    
    !Проинтегрируем диф. уравнение
    do while (z > dz)
        s = pi * (rshar**2 - abs(z - rshar)**2)
        dt = dz * s / (sotv * this%k * sqrt(2.0 * g * z))
        this%t = this%t + dt 
        z = z - dz
    end do  
 
  end subroutine method_t
end module telo


главная программа

Fortran
1
2
3
4
5
6
7
8
9
10
11
12
program fluid_outflow                           !Главная программа
use telo                                        !Подключим модуль
implicit none
 
    type (shar) :: object_1                     !Создадим объект шар
     object_1 % dotv = 0.12
     object_1 % vshar = 600
      
    call method_t (object_1)                    !Вычислим время
 
print*, "Время истечения жидкости = " , object_1 % t
end program fluid_outflow

Результат:
Код:
 Время истечения жидкости =    5924.95361

Однако можно решить задачу и другим, "более научным" способом, хотя сочетание "более научный" здесь конечно условно. Суть в том, что можно было не писать код решения диф. уравнения непосредственно в программе, а использовать для этого любую готовую библиотеку для численного решения диф. уравнений, к тому же предоставляющую какой-либо более точный метод решения. Такую библиотеку можно использовать универсально для решения разных задач. Главное - правильно ее применить. Имеем готовый модуль DU для решения дифференциального уравнения методом Рунге-Кутты (4-го порядка) - Милна.

модуль DU

Fortran
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
module DU
implicit none
 
public :: miln
private
 
 contains
 
!Программа численного решения задачи Коши для 
!дифференциального уравнения 1-го порядка методом Милна
subroutine miln (f, x1, x2, yfe, rez, eps, epsr, n, k)
implicit none
 
!Интерфейс к интегрируемой функции
interface
    real function f (x, y)
        real, intent (in) :: x, y
    end function f
end interface
 
!Описание переменных
real, dimension (0:5)   :: x, y, yk, z      !Рабочие массивы программы
integer, intent (out)   :: k                !Информация о достижении точности приближения Yk (0 или 1)
integer                 :: n, i             !n - Число шагов интегрирования, i - временная переменная 
real                    :: k1, k2, k3, k4   !Коэффициенты для метода Рунге-Кутты
real                    :: temp1, temp2     !Временные переменные
real                    :: h                !шаг интегрирования
real, intent (in)       :: x1               !Зачальное значение диапазона интегрирования
real, intent (inout)    :: x2               !Конечное знач. диапазона инт-я на входе, и на выходе суммир. по шагу  
real, intent (inout)    :: yfe              !начальное (на входе) и конечное (на выходе) значение Y
real, intent (out)      :: epsr             !Расчитанное максимальное отклонение Yk
real, intent (out)      :: rez              !Значение функции (производной)
real, intent (in)       :: eps              !Заданное значение точности приближения по скорректированному Y
 
!Решение
x = 0; y=0; yk=0; z = 0
    
    if (n < 50) n = 50
 
h = (x2 - x1) / float(n)
x(0) = x1
y(0) = yfe
 
    !Первые 4 шага вычисляются методом Рунге-Кутта 4-го порядка
    do i=1,4
        x(i) = x(i-1) + h
    end do
 
    yk(0) = y(0)
    z(0) = f(x(0), y(0))
    
    do i= 0, 3
        k1 = h * f(x(i), y(i))                  
        k2 = h * f(x(i) + h/2, y(i) + k1/2)     
        k3 = h * f(x(i) + h/2, y(i) + k2/2)
        k4 = h * f(x(i) + h, y(i) + k3)     
        y(i+1) = y(i) + (k1 + 2*k2 + 2*k3 + k4) / 6.0
        z(i+1) = f(x(i+1), y(i+1))
        yk(i+1) = y(i+1)        
    end do
        
    !Последующие шаги методом Милна  
    k = 0
    temp2 = 0.0
    x(5)  = x(4)
    temp1 = x(5)
        
    do i=5, n
        x(5) = temp1 + h
        temp1 = x(5)
        y(5) = yk(1) + 4*h/3 * (2*z(2) - z(3) + 2 * z(4))
        yk(5)  = yk(3) + h/3 * (z(3) + 4*z(4) + f(x(5), y(5))) 
        epsr = abs(y(5)-yk(5))/29.0
            if (epsr > eps) k = 1
            if (epsr > temp2) temp2 = epsr 
        z(5) = f(x(5), yk(5)) 
        x  = cshift (x,  1, 1) !Сдвиг массивов данных
        y  = cshift (y,  1, 1)
        yk = cshift (yk, 1, 1)
        z  = cshift (z,  1, 1)
    end do
        
!Результат вычислений
x2  = x(4) 
yfe = yk(4)
rez = z(4)
epsr = temp2
 
end subroutine miln
end module DU

Используя модуль остается написать код главной программы и саму функцию в соответствие с требуемым интерфейсом. Решаемая модулем задача - на задаваемом отрезке X[0...n] c начальным условием Y(0) при заданном числе шагов, и соответственно величине шага, производится интегрирование функции с поиском значения Y(n). Иначе говоря при определенных исходных данных численно находится значение неизвестной функции y=f1(x) в нужной нам точке, это и есть решение задачи Коши. Общий вид диф. уравнения, которое решает модуль:

http://www.cyberforum.ru/cgi-bin/latex.cgi?{y}^{'}=\frac{dy}{dx}=f(x, y)

При этом функция http://www.cyberforum.ru/cgi-bin/latex.cgi?f(x,y), записывается в виде отдельной программной единицы (либо размещается в коде главной программы, которую мы напишем). В нашем случае это правая часть решамого нами уравнения, приведенного в самом начале. Итак нам нужно написать функцию, задать начальные условия и "подсунуть" это модулю для решения. При этом важно правильно определить где X, а где Y. В нашем примере X - это высота поверхности жидкости, а точнее промежуток от максимума до минимума X[0...n] где за значение X(0) мы принимаем максимальный уровень жидкости, а за значение X(n) - пустой бак, Y - время, которое будет потрачено на вытекание жидкости, здесь точно так же фигурирует промежуток времени Y[0...n] при начальном значении Y(0)=0 и искомом значении Y(n).

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

главная программа

Fortran
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
program fluid_outflow
use DU
implicit none
real            :: x1, x2, y, rez, epsr
real, parameter :: eps = 1e-3
integer         :: n, k
 
!Пропишем интерфейс к интегрируемой функции
interface
    real function f (x, y)
        real, intent (in) :: x, y
    end function f
end interface
 
!Начальные параметры
x1 = 10.464 !Нулевая отметка (диаметр шара)
x2 = 1e-4   !Конечная отметка (уровень днища с отверствием)
y  = 0      !Значение времени в начальный момент
n  = 500    !Число шагов интегрирования
 
 !Решим диф. уравнение
 call miln (f, x1, x2, y, rez, eps, epsr, n, k)
 
!Выведем результат
print '(a, f9.2, a)', "t = ", y, " секунд"
end program fluid_outflow


Интерфейсный блок к функции в главной программе можно не прописывать, если код функции разместить в самой программе, отделив ее ключевым словом contains. Функция должна быть записана в соответствие с интерфейсом модуля (функция f и две входных переменных x и y, хотя переменные x и y в функции могут быть обозначены и другими буквами).

интегрируемая функция

Fortran
1
2
3
4
5
6
7
8
9
10
11
real function f (x, y)
implicit none
real, intent (in)   :: x, y
real                :: s
real, parameter     :: pi = 3.1416, g = 9.8067      !Константы
real, parameter     :: r = 5.232, so = 0.01130976   !Радиус шара и сечение отверстия
 
    s = pi * (r**2 - abs(x - r)**2)                 !Диф. уравнение (числитель, S сечения шара)
    f = - s / (so * sqrt(2.0 * g * x))              !Диф. уравнение (основное)
 
end function f

Результат:
Код:
t =  5924.13 секунд
Изменяя значения начального и конечного уровня жидкости можно расчитать другой промежуок времени, например если резервуар в исходной точке наполнен только наполовину.

Видим, что разница во времени, по сравнению с предыдущим методом расчета, получилась менее чем 1 секунда, при результате почти 6000 секунд такая разница не существенна. При этом нужно учесть, что в модуле применен более точный метод решения диф. уравнений, однако использовано более чем в 20 раз меньше шагов интегрирования - 500, вместо более чем 10 000, которые использовались для решения методом Эйлера.

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

PS.
Немного за рамками решения задачи остается вопрос поведения производной, но ведь как минимум интересно понять что она выражает и понаблюдать что с ней происходит по мере вытекания жидкости из резервуара. Как мы видим из первой формулы (левая ее часть) она выражает скорость изменения уровня жидкости в резервуаре. Модуль DU в своих выходных параметрах выдает значение производной, просто мы его не вывели, но если вывести несколько значений (обратных, мы ведь перевернули формулу) на разных шагах интегрирования, то мы увидим следующуую картину:

Код:
 
Производная     Высота

-0.236 		10.464 
-0.011 		10.004    
-0.00367	9.020 
-0.00230 	8.016  
-0.00174 	7.011    
-0.00146   	6.007    
-0.00131   	5.002    
-0.00124 	4.018    
-0.00123   	3.014    
-0.00133   	2.009    
-0.00168   	1.005    
-0.00226   	0.503 
-0.143     	0.0001
Отрицательный знак производной означает, что наша функция убывающая, ведь уровень жидкости в резервуаре уменьшается. Видно, что производная принимает самые маленькие по модулю значения на уровене, когда высота до поверхности жидкости составляет около 3-х метров (если точнее то около 3,5 м.), при диаметре резервуара 10.464 м. И это вполне объяснимо. Вспомним, что наш резервуар это шар, ближе к центру шара площадь его сечения максимальна, поэтому скорость изменения уровня жидкости, производная, минимальна. Минимальное значения не точно по центру, а именно ниже центра шара объясняются тем, что с понижением уровня жидкости, вследствие снижения напора, уменьшается и поток, представляющий собой произведение скорости истечения жидкости, определяемой формулой Торичелли, и сечения отверстия через которое истекает жидкость. Данное произведение записано в числителе правой части нашего диф. уравнения. В то же самое время максимальное по модулю значение производной можно увидеть в самом начале процесса, это то же хорошо укладывается в логику, ведь площадь поперечного сечения шара сверху невелика, а напор как раз максимален.

PS2
Интересно сравнить точность расчетов методами Милна и Эйлера на небольшом числе шагов. Видно, что график построенный по методу Милна идет очень ровно и при малых значениях шага дает значительно лучшее приближение. График построенный по методу Эйлера показывает меньшую точность. При совсем малом числе шагов, 10...20 (на графике не показано), величина ошибки расчета методом Эйлера достигает уже совершенно неприемлемых значений. Но провал в точности ожидаем, а вот неровность графика несколько озадачила, эта неровность присутствует и при более мелком шаге, хотя амплитуда "правалов и выбросов" становится много меньше.
Миниатюры
Нажмите на изображение для увеличения
Название: Miln-Eyler.png
Просмотров: 932
Размер:	19.7 Кб
ID:	5798  
Размещено в Без категории
Просмотров 248 Комментарии 0
Всего комментариев 0
Комментарии
 
КиберФорум - форум программистов, компьютерный форум, программирование
Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2020, vBulletin Solutions, Inc.