Добавил:
Upload Опубликованный материал нарушает ваши авторские права? Сообщите нам.
Вуз: Предмет: Файл:

caplin_nikulin_modelirovanie_v_metallurgii

.pdf
Скачиваний:
218
Добавлен:
13.03.2016
Размер:
10.25 Mб
Скачать

ва. Нарушение свойства транспортивности эквивалентно возникновению фиктивных (счетных) источников (стоков) в конечноразностном аналоге уравнения переноса.

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

5.4. Способы аппроксимации конвективных членов

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

Схема с центральными разностями

(uT )

ui +1Ti +1 ui 1Ti1

(5.47)

x

2hx

 

 

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

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

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

 

 

 

T

+1

T

 

 

 

ui

i

i

,

ui < 0,

 

 

 

 

 

(uT )

 

 

hx

 

 

 

Ti Ti1

 

 

(5.48)

x

 

 

 

ui

, ui 0.

 

 

 

 

 

 

 

hx

 

 

 

 

 

 

 

 

 

 

 

 

 

 

181

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

Вторая схема с разностями против потока известна как

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

(uT )

uпTп uлTл

,

(5.49)

x

 

 

hx

 

где uп = (ui +1 + ui )2, uл = (ui 1 + ui )2 , а значения Т выбираются в зависимости от знака усредненных скоростей:

Tп = Ti , uп 0,

Tл = Ti 1 , uл0,

Ti +1 , uп < 0.

Ti , uл < 0.

Схема с донорными ячейками обладает как свойством транспортивности, так и свойством консервативности. Формально она имеет первый порядок точности, однако усреднение скоростей сохраняет в ней кое-что от второго порядка точности. Поэтому схема (5.49) имеет меньшую по сравнению со схемой (5.48) величину счетной диффузии.

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

5.5. Аппроксимация граничных условий

Аппроксимацию граничных условий рассмотрим на примере граничных условий теплообмена 3-го рода для правой границы (рис. 5.8)

182

Рис. 5.8. Фрагмент сетки у правой границы

−λ

T

= α (T

T ) .

 

 

n

п

c

 

 

 

Применяя формулу односторонней разности, получим приближение:

−λ TN TN 1 = α (TT ) ,

hx

N C

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

 

T

+

 

λ

 

T

 

 

α

 

 

 

 

 

C

 

hx

 

 

TN = −

 

 

 

 

 

 

 

.

(5.59)

 

1 +

λ

 

 

 

 

 

 

 

α

hx

 

 

 

 

 

 

 

 

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

При работе со схемами второго порядка точности (Кранка – Николсона и др.) необходимо использовать более точную аппроксимацию граничных условий. Для этого запишем разложение температуры в ряд Тейлора в окрестности границы i = N:

T

= T

T

h

+

 

2T hx2

+ ... .

(5.60)

 

 

 

 

 

 

 

x2 2

N 1

N

 

x x

 

 

Учитывая три члена разложения, получим:

T

TN TN 1

+

2T hx2

+ ... .

 

 

 

 

 

 

 

x

hx

x2 2

 

183

Запишем вторую производную в конечных разностях:

T TN TN 1+ TN 2 2TN 1 + TN = TN 2 4TN 1 + 3TN .

x

hx

2hx2

2hx2

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

−λ

TN 2 4TN 1 + 3TN

= α

(T

T

)

 

 

2h2

N

C

 

 

x

 

 

 

получается более точное по сравнению с (5.59) значение температуры на границе:

 

T

+ (4T

 

T

)

λ

 

 

 

TN +1 =

c

N

 

N 1

2α hx2

 

 

 

 

 

 

 

(5.61)

 

1

+

 

3λ

 

 

 

 

 

2α hx2

 

 

 

 

 

 

 

 

 

 

 

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

5.6. Методы решения сеточных уравнений

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

ATi 1 + BTi + CTi+1 = Fi , i =1, 2, ..., N 1,

(5.62)

которую можно записать в векторно-матричном виде:

 

[H ] { T } = { F }

(5.63)

184

 

или

B

C

 

 

 

T

 

 

F

 

 

 

 

 

 

 

1

 

 

1

 

 

A B C

 

 

 

T2

 

 

F2

 

 

 

A B C

 

 

T3

 

F3

 

(5.64)

 

 

 

 

 

 

 

=

 

,

 

 

 

 

 

 

 

 

 

 

 

 

A

B

C

 

 

 

 

 

 

 

 

 

A

 

 

 

 

 

 

 

 

 

 

B

T

 

F

 

 

 

 

 

 

 

N 1

 

 

N 1

 

 

где [H]− матрица коэффициентов; {T} − вектор-столбец неизвестных значений искомого параметра Т в узловых точках; {F} − неизвестный вектор столбец, характеризующий краевые условия и распределение параметра Т на предыдущем временном слое.

Матрица [H] обладает рядом специальных свойств, которые необходимо использовать при решении системы уравнений. Она имеет высокий порядок, зависящий от густоты сетки, которая может достигать для современных компьютеров нескольких десятков тысяч. Матрица является редко заполненной с размещением ненулевых элементов по диагонали в три ряда. Такие матрицы называются ленточными трехдиагональными. Важным свойством является симметрия матрицы относительно ее диагонали, вытекающая из равенства коэффициентов А и С. Указанные свойства матрицы [H] позволяют занимать незначительное место для ее хранения в запоминающем устройстве компьютера, поэтому матрицу [H] называют порождающейся в отличие от хранящейся матрицы.

Перейдем к рассмотрению эффективных способов решения системы (5.62).

Метод прогонки

Метод прогонки является модификацией метода исключения Гаусса, учитывающей свойства матрицы H. Решение системы (5.64) в узловой точке ищется в виде линейной функции. В частности, для (i− 1)-ой точки эта функция имеет вид:

185

Ti 1 = β iTi + zi ,

(5.65)

где β i , zi − неизвестные пока вспомогательные коэффициен-

ты. Подставим (5.65) в (5.62):

A(β iTi

+ zi ) + BTi + CTi+1 = Fi ,

(5.66)

откуда находим:

 

 

 

 

 

 

 

T

= −

C

T

Azi Fi

.

(5.67)

 

 

i

 

Aβ i + B i+1

 

Aβ i + B

 

Полученное соотношение имеет ту же форму, что и функция (5.65), только для i-й точки

 

 

Ti

= β i+1Ti+1 + zi +1 ,

 

 

 

(5.68)

откуда заключаем, что

 

 

 

 

 

 

 

 

β

= −

 

C

;

z = −

Azi

Fi

.

(5.69)

 

 

 

 

 

i+1

Aβ

 

+ B

i +1

Aβ i

+ B

 

 

 

i

 

 

Полученные коэффициенты называются прогоночными коэффициентами, а формулы (5.68–5.69) дают процедуру решения.

Сначала при i = 1, 2,..., N– 1 считаются прогоночные коэффициенты (5.69), при этом начальные значения прогоночных коэффициентов β 1 , z1 определяются из граничных условий на левой границе. Эта операция называется прямой прогонкой. После определения всех β i , zi в обратном направлении (i = N− 1,..., 1) сучетом значения параметра TN , найденных из граничного условия на правой границе, по формуле (5.68) последовательно находятся неизвестные значения Ti в узловых точкахсетки.

Определение начальных значений прогоночных коэффициентов рассмотрим на примере граничных условий теплообмена. На левой границе условие конвективного теплообмена:

λ

T

= α (T

T )

(5.70)

 

 

x

п

c

 

 

 

 

 

186

может быть представлено в конечно-разностном виде (рис. 5.9):

λ

T1 T0

= α (T0 Tc ) .

(5.71)

 

 

h

 

Рис. 5.9. Фрагмент сетки у левой (а) и правой (б) границ

Отсюда находим:

 

 

 

λ

 

 

 

 

 

Tc

 

 

 

T0 =

 

 

α h

 

T1

+

 

 

.

(5.72)

 

 

 

 

 

 

λ

 

 

1 +

λ

 

 

 

1 +

 

 

 

 

α

h

 

 

h

 

 

 

 

 

 

 

α

 

 

Сравнивая эту формулу с соотношением (5.68) для левой

границы

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

T0 = β 1T1

+ z1 ,

 

 

 

(5.73)

получаем начальные значения прогоночных коэффициентов:

 

 

 

 

 

λ

 

 

 

 

 

 

 

Tc

 

β

1=

 

 

α

h

 

 

 

 

z1=

 

 

 

 

 

 

 

 

;

 

 

 

 

.

(5.74)

 

 

 

λ

 

 

 

 

λ

 

 

 

1 +

 

 

 

 

 

1 +

 

 

 

 

 

 

α

h

 

 

 

 

 

 

 

 

 

 

 

α h

 

Запишем условие теплообмена на правой границе:

 

 

−λ

T

= α

(T

T ) ,

(5.75)

 

 

 

 

 

 

 

 

x

п

c

 

 

 

 

 

 

 

 

 

 

 

 

 

 

в конечных разностях:

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

λ

TN TN 1

= α

(TN

Tc ) .

(5.76)

 

 

 

 

 

 

 

h

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

187

Отсюда находим:

 

λ

+1

 

α h

 

 

T =

α h

T .

(5.77)

 

 

T

 

 

 

λ

N

 

λ

 

 

N +1

 

c

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

α

h

 

 

 

 

Запишем соотношение (5.68) для правой границы:

 

TN 1 = β N TN + zN .

 

(5.78)

Приравнивая правые части (5.77), (5.78), получим искомое значение температуры на правой границе

 

 

 

 

 

 

 

 

 

 

 

 

 

 

λ

 

zN + Tc

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

TN =

 

 

бh

 

 

 

 

 

 

.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

1 +

λ

(1 вN )

 

 

 

 

 

 

 

 

 

 

 

 

 

бh

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Запишем алгоритм метода прогонки:

 

 

 

 

 

λ

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Tc

 

 

 

 

 

β

1=

 

 

α h

 

 

 

 

 

z1=

 

 

 

 

 

 

 

 

 

;

 

 

 

;

 

 

 

 

 

 

 

λ

 

 

 

 

 

λ

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

1 + α h

 

 

 

 

 

 

1 +α

h

 

 

 

 

 

 

 

 

 

 

 

 

 

C

 

 

 

 

 

 

 

 

 

 

 

Azi Fi

 

β

= −

 

 

 

 

 

 

 

;

 

 

 

z

= −

 

 

 

,

 

 

 

 

 

 

 

 

 

 

 

 

 

 

i+1

 

 

Aβ i+ B

 

 

 

 

i+1

 

 

Aβ +i B

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

i = 2,

3,

 

...,

N;

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

λ

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

α h zN + Tc

 

 

 

 

 

 

 

 

 

 

TN =

 

 

 

 

 

;

 

 

 

 

 

 

 

 

 

 

 

λ

 

 

(1

− β N )

 

 

 

 

 

 

 

 

1 +

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

α

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

h

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Ti = β i+1Ti +1+ zi +1 ,

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

i = N 1,

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

...,

0.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

(5.79)

(5.80)

Пример программы, реализующий алгоритм метода прогонки, приведён в лабораторной работе № 3.

188

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

Метод последовательной линейной верхней релаксации

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

Рассмотрим один из эффективных итерационных методов – метод последовательной линейной верхней релаксации, итерационная процедура которого применительно к разностному уравнению (5.62) имеет вид

Ti q =

γ

(Fi ATi q1 CTi+q11 ) + (1− ω )Ti q1 , i= 1, 2, , N1, (5.82)

B

 

 

где q

номер итерации; γ – параметр релаксации. При γ =1

получаем процесс последовательных смещений, или процесс Зейделя. Введение параметра верхней релаксации 1 ≤ γ ≤ 2 по-

зволяет ускорить сходимость итерационного процесса (21), причем наибольшая скорость сходимости имеет место при оптимальном значении параметра релаксации γ = γ опт . Последнее

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

γ опт =

 

 

 

2

 

 

 

.

(5.83)

 

 

π

 

 

π

 

1 +

 

sin

 

 

sin

 

2

 

 

 

 

 

 

 

 

 

2N

 

2N

 

 

 

 

 

 

 

 

 

 

189

Эта формула применима для одномерной области с регулярной сеткой. Расчет по формуле (5.82) с учетом (5.83) продолжается до тех пор, пока искомое решение не будет удовлетворять требуемой точности:

1

T q1

 

 

≤ ε ,

 

i

 

 

(5.84)

T q

 

 

 

 

 

 

i

 

 

max

 

 

 

 

где ε – требуемая точность.

5.7. Алгоритм решения сопряженных уравнений теплообмена

Рассмотрим в общих чертах процедуру решения сопряженной задачи конвективного теплообмена, в которой совместно решаются уравнения конвекции, например, в переменных завихренность и функция тока (ω ψ ) и уравнение переноса энергии для температуры Т. Рассмотрим вычислительный цикл для нестационарных уравнений ω ψ T-системы (рис. 5.10).

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

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

190

Соседние файлы в предмете [НЕСОРТИРОВАННОЕ]