Метод контрольных объемов (1D)

Метод контрольных объемов часто используется при решении задач, связанных с вычислительной гидродинамикой и теплообменом. Рассмотрим одномерный случай установившегося распределения температуры $T$ в теплоизолированном по длине стержне, на концах которого поддерживается постоянная температура.

$A$ — площадь поперечного сечения равна 0,001 кв.м, $K$ — коэффициент теплопроводности равен 1000 Вт/мК — это довольно много, коэффициент теплопроводности, например, алюминия равен 210-220 Вт/мК. Буквой P обозначается узел, для которого формируется уравнение, большими буквами W (west) и E (east) обозначаются западный и восточный соседние узлы. Маленькими буквами w и e обозначаются границы контрольного объема с узловой точкой P.

Метод контрольных объемов (1D) — общие сведения

Обозначения приняты в соответствии с книгой H.K. Versteeg, W. Malalasekera — An Introduction to Computational Fluid Dynamics. The Finite Volume Method (2-nd ed) см. стр.115-121.

Есть ещё две статьи по одномерному МКО. Перевод разделов 4.1, 4.2 здесь, а перевод раздела 4.3 с реализованным в данной статье примером здесь.

Реализация МКО в Matlab

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

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

$\begin{pmatrix} -300 & 100 & 0 & 0 & 0 \\ 100 & -200 & 100 & 0 & 0 \\ 0 & 100 & -200 & 100 & 0 \\ 0 & 0 & 100 & -200 & 100 \\ 0 & 0 & 0 & 100 & -300 \end{pmatrix} \cdot \begin{pmatrix} T_{1} \\ T_{2} \\ T_{3} \\ T_{4} \\ T_{5} \end{pmatrix} = \begin{pmatrix} -2 \times10^{4} \\ 0 \\ 0 \\ 0 \\ -1 \times10^{5} \end{pmatrix}$

В реализации не учитывается предположение о равномерности сетки, поэтому введены переменные $\delta_{w}$ и $\delta_{e}$ (по коду deltaXw и deltaXw), которые хоть и не отличаются в рассматриваемой задаче, могут быть различными, если сетка не равномерная.

Учет температурных ГУ сводится к тому, что значение температуры берется не в узловой точке (как в 4.9), а в граничной. Расстояние тоже берется до граничной точки (считаем, что узлы в середине контрольных объемов). Первый и пятый контрольные объемы являются граничными и обрабатываются отдельно от внутренних. Для пятого узла, например, с учетом неравномерности сетки, вместо (4.19) имеем:

$k A \left( \frac{T_{B}-T_{P}}{dx_{e}/2} \right) - k A \left( \frac{T_{P}-T_{W}}{dx_{w}}\right)=0$

Здесь есть узловые точки контрольных объемов и граничные точки. Отличием МКО от метода конечных элементов является то, что ГУ можно задать не в узловой точке контрольного объема, а в граничной точке.

Уравнения системы записываются для внутренних и для граничных контрольных объемов. 5 узловых точек, значит 5 контрольных объемов (три внутренних и два граничных), и уравнений в системе тоже 5. Однако точек в сетке — 11. Это с учетом граничных точек элементов. Поэтому реализация содержит два массива npts (узлы контрольных объемов) и bpts (границы контрольных объемов).

Код (MATLAB):


function RunCVM()
% Стр.118 Versteeg, Malalasekera - Introduction to CFD. The finite
% volume method. 2-nd ed.
% Коэффициент теплопроводности.
k = 1000;
% Площадь сечения стержня.
A = 10e-3;
% Температуры на краях стержня.
TA = 100;
TB=500;
% Кол-во контрольных объемов.
ncv = 5;
% Точки границ объемов.
bpts = (0:(0.5/ncv):0.5)';
% Создание массива узловых точек контрольных объемов.
npts = zeros(ncv,1);
% Инициализация массива узловых точек контрольных объемов.
for i = 1:ncv
npts(i) = bpts(i) + ((bpts(i+1)-bpts(i))/2);
end
% Печать созданной сетки.
plot(npts,1, 'bo',...
'MarkerSize',8,'MarkerFaceColor', 'w',...
'LineWidth', 2.0);
grid on;
hold on;
plot(bpts,1, 'rx',...
'MarkerSize',8,'MarkerFaceColor', 'w',...
'LineWidth', 1.0);
% Создание матрицы системы дискретного аналога.
a = zeros(ncv,ncv);
% Создание вектора правой части.
b = zeros(ncv,1);

% Расчет коэффициентов для каждого внутреннего контрольного объема.
for i=2:ncv-1
deltaXe = npts(i)-npts(i-1);
deltaXw = npts(i+1)-npts(i);
ae = (k/deltaXe)*A;
aw = (k/deltaXw)*A;
ap = -(ae+aw);
a(i,i) = ap;
a(i,i-1) = ae;
a(i,i+1) = aw;
end

% Коэффициенты для первого контрольного объема (с ГУ). Нет aw.
deltaXw = (npts(2)-npts(1))/2; % КО на границе, поэтому /2.
deltaXe = npts(3)-npts(2);
ae = k*A/deltaXe;
Sp = -k*A/deltaXw;
Su = k*A*TA/deltaXw;
ap = -(ae-Sp);
a(1,1) = ap;
a(1,2) = ae;
b(1) = -Su;

% Коэффициенты для последнего контрольного объема (с ГУ). Нет ae.
deltaXw = npts(2)-npts(1);
deltaXe = (npts(3)-npts(2))/2;
aw = k*A/deltaXw;
Sp = -k*A/deltaXe;
Su = k*A*TB/deltaXe;
ap = -(aw-Sp);
a(ncv,ncv) = ap;
a(ncv,ncv-1) = aw;
b(ncv) = -Su;

% Решение системы.
T = a\b;
% Вывод результата на экран.
fprintf("TA=%d\n\n",TA);
disp(T);
fprintf("TB=%d\n",TB);
end

Результаты

Метод контрольных объемов (1D) - Результат запуска программы. 5 контрольных объемов. — Метод контрольных объемов (1D) — Результат запуска программы. 5 контрольных объемов. Узлы контрольных объемов показаны синими кружочками, их границы показаны красными крестиками.

Метод контрольных объемов (1D) — общие сведения

Реализация МКО в Matlab

Результаты

Добавить комментарий Отменить ответ