Барилович В.А. Основы термогазодинамики двухфазных потоков и их численное решение

110

(

)

[

]

2

1

жжж1п

4

d

ssd

dx

dp

ww

π

τ

π

τ

+

−

−=

. (*)

Степень сухости парожидкостного потока определим по уравнению

dx

Qd

rG

dx

v

pc

см

1

&

=

.

Скорость жидкости находится из уравнения сплошности

(

)

жж

см

ж

1

f

Gx

с

v

ρ

&

−

=

.

Уравнения для греющего потока записываются так же, как в гомогенной

модели.

Предложенная физическая модель позволяет использовать традиционные

критериальные уравнения теплоотдачи и трения, полученные для однофазных

сред, и отказаться от уравнений вида [57]

⎥

⎦

⎤

⎢

⎣

⎡

⎟

⎠

⎞

⎜

⎝

⎛

⎟

⎠

⎞

⎜

⎝

⎛

⎟

⎠

⎞

⎜

⎝

⎛

⎟

⎠

⎞

⎜

⎝

⎛

−

⎟

⎠

⎞

⎜

⎝

⎛

−

+=

qР

m

F

x

CF

ж

г

ж

г

p

pжтр

Pr

1

μ

αα

,

которые справедливы лишь при определенных режимных параметрах.

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

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

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

скажется на теплообмене. Рост скорости пара усилит этот эффект. Таким

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

позволит

определить максимальную длину парогенерирующего канала.

Если учитывать скольжение между фазами, то система уравнений примет

вид:

⎥

⎦

⎤

⎢

⎣

⎡

−−−−−=

−

dx

Gd

wwbs

dx

dp

f

G

dx

dw

w

пж

жп1пжппп

п

)(

1

&

ττ

,

где

dx

Qd

rdx

Gd

рс

пж

1

&

=

−

,

⎥

⎦

⎤

⎢

⎣

⎡

+−−=

1пжжжж

ж

1

bs

dx

dp

f

G

dx

dw

w

ττ

&

,

dx

Gd

dx

dw

f

dx

df

w

пж

ж

1

−

−=+

&

ρ

,

.0

п

пп

п

пп

п

пп

ж

жж

ж

жж

=++++

dx

df

w

dx

dp

d

wf

dx

dw

f

dx

df

w

dx

dw

f

ρ

ρρρ

Разрешая ее относительно dp/dx , получим

111

1

п

2

ж

п

3

ж

п

ж

п

2

п

×

⎥

⎦

⎤

⎢

⎣

⎡

−

⎟

⎠

⎞

⎜

⎝

⎛

⎟

⎠

⎞

⎜

⎝

⎛

−

=

f

ρ

w

G

/dp)(dρ

w

dx

dp

п

&

⎥

⎦

⎤

⎢

⎣

⎡

⎟

⎠

⎞

⎜

⎝

⎛

−−

⎟

⎠

⎞

⎜

⎝

⎛

−−+×

−

dx

Gd

w

wρ

bτ

w

wρ

w

sτ

wρ

w

пж

п

ж

жж

пп

1пж

п

2

жж

пп

п

пп

жж

2

жж

пп

1

&

Зная dp/dx, найдем:

dxdw /

п

,

dxdw /

ж

,

dx

Gd

wdxw

dwf

dx

df

пж

жжж

жжп

1

−

+=

&

ρ

,

пж

fFf −=

.

Далее определим

5,0

1111

))((2 hhdb −=

,

1

11

1

п

24

b

d

hb

d

f

s −+=

и

п1ж

sds −=

π

.

Результаты расчета парогенерирующего канала представлены на рис.7.6.

Система уравнений была реализована в среде Visual Basic.

Рис. 7.6. Расчетное изменение температуры и скорости азота и воздуха в парогенерирующем канале.

Рисунок иллюстрирует ход кривых скорости и температуры агентов вдоль

канала на экономайзерном, парогенерирующем и участке перегрева пара. На

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

112

Раздел II. Движение частиц и двухфазного потока

в криволинейных каналах

В разделе I рассматривалось движение двухфазного и вскипающего

потоков в каналах переменного сечения с прямой осью. Однако во многих

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

однокомпонентных или многокомпонентных гетерогенных потоков в

криволинейных каналах (решетки турбомашин, колена трубопроводов и

т.д.) с

частичным или полным инерционным выпадением отдельных компонент на

стенках канала. В связи с вышесказанным возникает необходимость в изучении

таких потоков.

8. Движение капли в сепараторе влаги

Пар, движущийся в паропроводах энергетических систем, когда его

параметры близки к состоянию насыщения, из-за теплоотдачи в окружающую

среду становится влажным

. Однако, например, работа турбины на влажном

паре снижает ее внутренний относительный к.п.д., приводит к эрозионному

износу лопаток, уменьшает срок ее эксплуатации. Для борьбы с этими

негативными явлениями применяют сепараторы влаги. Существует множество

типов сепараторов влаги и частиц, но большинство из них основаны на

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

коаксиальными

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

результате возникающих при движении по окружности ускорений, частицы

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

пленку, которая стекает в нижнюю часть сепаратора. Отсепарированная влага

удаляется через патрубок. Для предотвращения срыва жидкости с поверхности

пленки паровым потоком (т.е. вторичного

увлажнения потока) скорость

входящего в сепаратор пара не должна быть слишком высокой, обычно она

составляет

с

Пϕ 0

≈20...30 м/с. С этой же целью в нижней части сепаратора

устанавливается коническая крышка. На рис.8.1 показана возможная

конструкция противоточного сепаратора влаги. Чтобы понять принцип

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

материальной точки (см. рис.8.2).

В полярной системе координат абсолютная скорость точки

М, движущейся в

плоскости по траектории

АВ, может быть разложена на относительную

скорость

•

==

ρ

τ

ρ

d

с

, направленную вдоль радиус-вектора

ρ

, и на переносную

•

==

ϕρ

τ

ϕ

ρ

ϕ

d

с

, перпендикулярную радиус-вектору. Так как скорости с

ρ

и с

ϕ

взаимно перпендикулярны, то

22

ϕρ

ссс +=

, а

() ()

c

cc

c

cc

ϕ

ρ

==

∧∧

cos,cos (см.

рис. 8.2). Абсолютное ускорение точки складывается из относительного

113

ускорения

••

==

ρ

τ

ρ

d

dc

a

, направленного по радиусу-вектору, переносного

e

a

r

,

которое состоит из переносного касательного

τ

e

a

r

и переносного нормального

en

a

r

,

enee

aaa

r

+=

τ

,

2

•••

+=

⎟

⎠

⎞

⎜

⎝

⎛

+=

ϕρϕρ

τ

ϕ

ρ

τ

ϕ

ρ

d

a

e

, и кориолисова ускорения

ca

c

r

×=

ϖ

2 , которое по направлению совпадает с

τ

e

a

r

,

••∧

=

⎟

⎠

⎞

⎜

⎝

⎛

=

ρϕϖ

τ

ϕ

ρρ

2sin2 cc

d

a

c

rr

.

Абсолютное ускорение точки

складывается из относительного

ускорения

••

==

ρ

τ

ρ

d

dc

a

, направленного по

радиусу-вектору, переносного

e

a

r

,

которое состоит из переносного

касательного

τ

e

a

r

и переносного

нормального

en

a

r

,

enee

aaa

r

+=

τ

,

2

•••

+=

⎟

⎠

⎞

⎜

⎝

⎛

+=

ϕρϕρ

τ

ϕ

ρ

τ

ϕ

ρ

d

a

e

, и кориолисова

ускорения

ca

c

r

×=

ϖ

2 , которое по направлению совпадает с

τ

e

a

r

,

••∧

=

⎟

⎠

⎞

⎜

⎝

⎛

=

ρϕϖ

τ

ϕ

ρρ

2sin2 cc

d

a

c

rr

. Из рисунка 8.3 видно, что проекция ускорения

точки на радиальное направление равна

ρτ

ρϖ

τ

ϕρρ

ρρ

2

c

d

dc

d

dc

aaa

en

−=−=−

′

= .

На направлении, перпендикулярном радиус-вектору, будем иметь

=+

⎟

⎠

⎞

⎜

⎝

⎛

=+=+=+=

ρρτ

ρ

ρτ

ϖ

ρϖ

τ

ϖ

ρ

ρϕϕρϕ

ρτϕ

ccc

d

cc

d

c

d

aaa

ce

222

.22

2

ρτρρτρρ

τ

ρ

τ

ρρ

ρϕϕρϕρϕϕρϕ

ϕ

cc

d

dccccc

d

dccc

d

c

d

dc

+=+−=+

⎟

⎠

⎞

⎜

⎝

⎛

−

=

Рис. 8.1

Рис. 8.2

114

Таким образом, мы получили, что

ρτ

ϕρ

ρ

2

c

d

dc

a −=

,

ρτ

ρϕϕ

ϕ

cc

d

dc

a +=

.

В цилиндрической системе координат уравнения

примут вид (см. рис.8.4)

,,,

•••

=== zccc

z

ϕρρ

ϕρ

222

z

cccc ++=

ϕρ

,

() ()

()

,cos

,cos,cos

c

z

zc

c

•

∧

•

∧

•

∧

=

==

ρ

ϕρ

ϕ

ρ

,,,2,

222

2

zz

aaaazaaa ++==+=−=

•••••••••

ϕρϕρ

ϕρϕρϕρρ

() () ()

a

z

za

a

••

∧

••••

∧

•••

∧

=

+

=

−

= cos,

2

cos,cos

2

ϕρϕρ

ϕ

ϕρρ

ρ

.

Перейдем к рассмотрению движения капли в сепараторе влаги. Согласно

закону Ньютона, сила инерции уравновешивается силами, приложенными к

телу

∑

=

n

i

FaG

1

r

. Так как сила аэродинамического сопротивления намного

больше массовых сил в поле тяжести и сил давления, то в координатной форме

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

()

,

8

отн

2

ccc

D

c

r

c

d

dc

G

rrx

u

r

Δ−=

⎟

⎠

⎞

⎜

⎝

⎛

−

ΚΠΠ

Κ

ρ

π

τ

()

,

8

отн

2

ccc

D

c

r

cc

d

dc

G

uux

ruu

Δ−=

⎟

⎠

⎞

⎜

⎝

⎛

+

ΚΠΠ

Κ

ΚΚΚ

Κ

ρ

π

τ

()

,

8

отн

2

ccc

D

c

d

dc

G

zzx

z

Δ−=

ΚΠΠ

ΚΚ

Κ

ρ

π

τ

где

6

3

Κ

ΚΚ

=

D

G

π

ρ

- масса капли;

()()()

222

отн zzuurr

ccccccc

ΚΠΚΠΚΠ

−+−+−=Δ -

относительная скорость.

С учетом уравнения линии тока

zur

c

dz

c

rd

c

dr

ΚΚΚ

==

ϕ

преобразуем уравнения

движения к виду

Рис. 8.3

Рис. 8.4

115

()

.

4

3

,

4

3

,

4

3

отн

2

отн

u

zzx

z

r

u

uux

u

r

u

r

rrx

r

c

cc

D

c

dz

dc

c

cc

D

r

c

d

dc

rc

c

cc

D

c

dr

dc

Κ

ΚΠ

ΚΚ

ΠΚ

Κ

ΚΠ

ΚΚ

Π

Κ

ΚΠ

ΚΚ

ΠΚ

Δ

−=

−

Δ

−=

+

Δ

−=

ρ

ϕ

ρ

Полагая

с

Пr

=0, получим

ΚΚ

Π

Κ

Δ

−=

D

c

rc

c

dr

dc

x

r

u

r

ρ

отн

2

4

3

. (*)

При входе влажного пара в сепаратор

с

Кr

=0, следовательно, ∞=

Κ

dr

dc

r

и

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

реализации. Однако если записать их в виде производных по

ϕ

, проблема

исчезает. Действительно,

.1

4

3

,1

4

3

,

4

3

отн

u

z

x

z

r

u

x

u

r

xu

r

c

cc

c

D

r

c

d

dc

cc

c

D

r

c

d

dc

c

cc

D

r

cc

d

dc

Κ

Π

ΚΚ

ΠΚ

Κ

Π

ΚΚ

Π

Κ

ΚΚ

Π

Κ

Δ

⎟

⎠

⎞

⎜

⎝

⎛

−=

−Δ

⎟

⎠

⎞

⎜

⎝

⎛

−=

Δ

−=

ρ

ϕ

ρ

ϕ

ρ

ϕ

При расчете движения

капель в сепараторе считалось,

что ,/

0

rrcc

cp

⋅=

ΠΠ

ϕϕ

где

r

cp

=(r

нар

+r

вн

)/2, а с

Пr

=0. Влияние

капель на паровой поток не

учитывалось,

с

Пz

=const.

На основании приведенных

формул составлена программа

SEPARATOR, позволяющая

рассчитывать движение капли в

сепараторе. Результаты расчета

представлены на рис.8.5, текст программы - в приложении №6.

9. Движение двухфазного потока в криволинейном канале

В предыдущем параграфе рассмотрено движение капли в сепараторе

влаги, определялась точка встречи капли со стенкой, движение же пленки,

образовавшейся в результате выпадения капель на стенку, не рассматривалось.

Решение задачи о течении парокапельного потока и пленки в каналах в

трехмерной постановке представляет сложную научную проблему

диссертационного плана, поэтому ограничимся рассмотрением движения

Рис. 8.5

116

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

межлопаточного канала гидропаровой турбины (ГПТ), работающей на

вскипающем потоке. Использование ГПТ в геотермальной энергетике вместо

расширителей и в энергетических системах с глубокой утилизацией теплоты

повышает тепловую эффективность установок.

Будем считать, что на входе в канал высоковлажный парокапельный поток

имеет неравновесную структуру (

Т

к

>Т

п

,

пк

ww

r

≠

), несущий паровой поток

полностью заполняет канал, радиальная составляющая скорости парового

потока

w

пr

=0, а тангенциальная проекция зависит от радиуса по закону

const

п

=⋅ rw

ϕ

, где r- текущий радиус канала, срыв жидкости с поверхности

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

влияние на поток касательного напряжения трения на стенке

τ

w

и силу

сопротивления капель, которая снижает среднюю скорость несущего парового

потока на 15…17%. Для того чтобы ответить на вопрос, правомерно ли

использование модели плоского канала при наличии вращения, когда на поток

действуют центробежные и кориолисовы силы, необходимо оценить величину

Δ

h перемещения жидкой частицы по корыту пера лопатки в направлении оси h

(см. рис. 9.1).

Воспользуемся уравнением (см.

§10, ф-ла 10.7а)

∑

+++=

hi

u

F

h

w

wh

d

hd

,

2

пл.

2

ωω

τ

, (9.1)

где

ω

- угловая скорость рабочего колеса, рад/с; h – расстояние от оси

вращения, м;

u

w

пл.

- проекция скорости жидкостного элемента на ось и.

Пренебрегая всеми внешними силами

hi

F

,

в направлении h, кроме

центробежных и кориолисовых сил, оценим величину Δ

h за время Δ

τ

нахождения жидкой частицы в рабочем канале на основании опытных данных

одного из режимов работы ГПТ-100 мощностью 100 кВт:

u

w

пл.

=33 м/с,

ω

=397 с

-1

, l

к

=26,3 мм (длина корыта лопатки вдоль образующей), при этом

3

пл.к

10797,0/

−

⋅==Δ

u

wl

τ

с.

Записав уравнение второго порядка (9.1) в виде системы двух обыкновенных

дифференциальных уравнений первого порядка

h

w

wh

d

dw

w

d

dh

u

h

2

пл.

2

2и ++==

ωω

ττ

,

где

w

h

– проекция скорости жидкостного элемента на ось h, и решив систему

численным методом (см. программу RK1_GPT), находим, что Δ

h=3,03 мм. При

расчете учитывалось, что во время движения частиц по корыту лопатки

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

2

τ

Δ

.

В действительности, величина и направление кориолисова ускорения меняются

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

117

жидкой частицы

пл

w

r

, а, следовательно, меняется кориолисово ускорение

()

пл

wa

r

×=

ω

2

к

. Проведенная оценка Δh для данных условий позволяет с

достаточной точностью проводить расчет движения потока в межлопаточном

канале на основе модели плоского канала.

u

ϕ

d

ϕ

u

z

h

w

r

Рис. 9.1

Получим уравнения сохранения для стационарного режима течения в

полярной системе координат [50-52].

Изменение секундной массы в фиксированном паровом элементе

обусловлено потоками массы от испаряющихся капель (

Т

к

>Т

п

), тогда

(

)

dVnjDfwd

vп-к

2

кппп

πρ

ϕ

=

, (9.2)

где

dV – объем ячейки. Так как

ϕ

Frdd

V

=

(F – площадь поперечного сечения

ячейки,

d

ϕ

- элементарный угол поворота потока), то (9.2) можно переписать в

виде

()

rnjDw

d

vп-к

2

кппп

πϕρ

ϕ

= , (9.3)

где

Ff

пп

=

ϕ

- истинное паросодержание. После дифференцирования

уравнения сплошности (9.3) имеем

ϕ

ϕϕ

ρ

ρϕρ

π

ϕ

ϕϕϕ

d

w

d

dp

d

w

rnjD

d

dw

v

п

пп

п-к

2

к

п

−−= . (9.4)

Чтобы исключить в уравнении (9.4) производную

ϕ

dd

п

, найдем связь между

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

1 или,

кп

=++=

dV

dVdVdV

, вводя обозначения

6

3

к

v

nD

dV

π

γ

== и

dV

п

=

β

, получим

6

1

3

к

п

v

nD

π

β

−= .

Тогда

118

ϕ

π

ϕ

β

d

dD

n

D

d

v

к

2

кп

2

−=

. (9.5)

Расход капельного потока

6

3

к

кк

D

nG

π

ρ

&

= ( n

&

- поток капель, 1/с) после

дифференцирования по углу можно записать в виде

.

22

к

2

к

кк

к

2

к

ϕ

π

ρ

ϕ

π

ρ

ϕ

d

dDD

Fwn

d

dDD

n

d

Gd

v

==

&

(а)

С другой стороны,

rFnjD

d

Gd

vп-к

2

к

π

ϕ

−=

&

. (b)

Приравнивая (а) и (b), определим производную

r

w

j

d

dD

кк

п-кк

2

ρϕ

ϕ

−= . (9.6)

Перепишем (9.5) с учетом (9.6)

кк

п-к

2

к

п

ρ

π

ϕ

β

ϕ

w

rjnD

d

v

=

. (9.7)

Можно показать, что

()

(

)

2/

ппп

ϕϕϕ

ϕ

β

d+

+

= , а 2/

2

ппп

ϕϕ

ϕϕβ

ddd +=

или

ϕ

β

пп

dd ≈ . Тогда уравнение сплошности примет вид

ϕ

ρ

ρρρϕ

π

ϕ

ϕϕ

d

dp

d

w

rjnD

d

dw

v

п

кк

п

пп

п-к

2

к

п

1

−

⎟

⎠

⎞

⎜

⎝

⎛

−= . (9.8)

Для нашего случая

0,0,0

п

===

ϕτ

ϕ

d

dr

w

d

dw

r

и

rr

pp

dr

d

r

p

d

dp

⎟

⎠

⎞

⎜

⎝

⎛

∂

=

⎟

⎠

⎞

⎜

⎝

⎛

∂

+

⎟

⎠

⎞

⎜

⎝

⎛

∂

=

ϕϕ

ϕ

.

Уравнение движения парового потока запишется так:

пп

п-к

2

к

п

к

п

к

п

2

к

пп

п

11

8

1

ϕρ

π

ϕ

π

ϕρϕ

ϕ

rjnD

w

wrn

w

D

с

d

dp

wd

dw

v

vx

⎟

⎠

⎞

⎜

⎝

⎛

−−Δ

⎟

⎠

⎞

⎜

⎝

⎛

−−−=

, (9.9)

где

()

2

кп

2

кп

ϕϕ

wwwww

rr

−+−=Δ .

Приравнивая правые части (9.8) и (9.9), определим градиент давления по углу

вдоль канала

.11

1

8

1

пп

п-к

2

к

кк

пп

п

к

п

к

п

2

к

п

2

п

пп

⎭

⎬

⎫

⎥

⎦

⎤

⎢

⎣

⎡

⎟

⎠

⎞

⎜

⎝

⎛

−+

⎟

⎠

⎞

⎜

⎝

⎛

−+

+

⎪

⎩

⎪

⎨

⎧

Δ

⎟

⎠

⎞

⎜

⎝

⎛

−

⎟

⎠

⎞

⎜

⎝

⎛

−

=

ϕρ

π

ρ

ϕ

π

ρ

ϕ

rjnD

w

wrn

w

D

с

ddp

w

d

dp

v

vx

(9.10)

Для испаряющейся капли уравнение движения в координатной форме с

учетом уравнения линии тока

ϕ

кк

w

rd

w

dr

r

=

при 0

п

=

r

w имеет вид

119

⎪

⎩

⎪

⎨

⎧

−=

−

⎟

⎠

⎞

⎜

⎝

⎛

−=

,

~

,1

~

к

пк

r

w

Fw

d

dw

wr

w

F

d

dw

rr

r

ϕ

ϕϕ

ϕ

(9.11)

где

кк

п

4

3

~

D

w

cF

x

ρ

Δ

=

;

()

к

Refc

x

= ;

п

пк

к

Re

μ

ρ

wD

Δ

=

.

Разрешая уравнение количества движения испаряющегося пристенного слоя

жидкости переменной массы (см. формулы 3.37, 3.43) относительно

ϕ

dwd

пл

,

получим

()

плплплплрлпл-прлпл

oc

рл

плк

пл

)(

ρδδττ

ϕ

δ

ϕϕ

ϕ

ϕϕ

wRR

d

dp

d

Gd

h

ww

d

wd

w

⎟

⎠

⎞

⎜

⎝

⎛

−+−−

−

=

&

, (9.12)

где

h

рл

– высота канала, R

рл

– радиус корыта рабочей лопатки.

При расчете пристенного слоя жидкости капля считалась выпавшей, если

ее траектория пересекает поверхность слоя, при этом слой получает секундную

массу жидкости

6

к

3

кoc

ϕρπ

dDnGd

′

=

&

, где

яч

znn

&&

=

′

- элементарный поток капель,

проходящий через одну ячейку первого по углу ряда. Начало траекторий

капель, имеющих размер

D

к

и условно содержащих массу всех пересекающих

ячейку капель – в серединах граней ячеек на входе в расчетную область.

Из уравнения сплошности испаряющегося пристенного слоя при

отсутствии срыва жидкости с его поверхности находим

()

ϕ

δ

ρδ

ϕϕ

δ

d

wd

w

wRj

d

Gd

hd

d

пл

плплплрлп-пл

oc

рл

пл

1

−

⎥

⎦

⎤

⎢

⎣

⎡

−−=

&

. (9.13)

Записав уравнение теплового баланса между каплей и паровым потоком,

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

()()

п-кпкп-к

кккк

к

~

6

jrTT

wcD

r

d

dT

p

+−−=

α

ρϕ

ϕ

, (9.14)

где

r

~

- удельная теплота парообразования, Дж/кг. Производную

ϕ

d

dT

пл

найдем

из уравнения энергии с учетом материального баланса, составленного для

фиксированного элемента пристенного слоя

()

[]

()}

.

1

пл

рлплрл

пплп-плп-пл

*

плпл

*

п

oc

*

пл

*

к

плпл

пл

ϕ

δ

α

ϕϕ

d

wd

c

w

hR

TTjiTi

d

Gd

ii

Gc

d

dT

p

−−×

×

⎪

⎩

⎪

⎨

⎧

−+−

′′

−−=

&

(9.15)

Концентрация капель в ячейке определяется выражением

рл

1

к

hdrdr

w

dln

n

i

ii

v

ϕ

∑

=

′

=

&

, (9.16)

где

dl

i

– длина траектории единичной капли внутри ячейки.

Барилович В.А. Основы термогазодинамики двухфазных потоков и их численное решение

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