Институт механики им. Р.Р.Мавлютова УНЦ РАН

Важным звеном топливно-энергетического комплекса России являются магистральные нефтепроводы (МН). Одной из сложных задач в процессе эксплуатации МН является определение оптимальных режимов работы насосных агрегатов, с точки зрения минимума затрат электроэнергии в фактическом или денежном выражении. Сложность задачи заключается в большом объеме возможных режимов работ, зависящих от многих параметров, таких как: физико-химические свойства нефти, схемы работ, состояние и количество находящегося в работе насосно-силового оборудования. Так для примера, для МН работающим по схеме «из насоса в насос» и включающего в технологический участок («от емкости до емкости») шесть нефтеперекачивающих станций (НПС) количество возможных режимов работы составляет несколько сот тысяч. Определить из этого числа ряд оптимальных довольно сложная задача. Проблема усугубляется технологическими особенностями процесса эксплуатации МН. К ним относятся: возможные ограничения давлений по трассе МН, возможность работы на разных технологических схемах, временные ограничения в работе насосно-силового оборудования находящегося в ремонте, ограничения по объемам перекачки, изменение допустимых давлений на НПС, большая протяженность МН и т.д. Например ОАО «АК «Транснефть» эксплуатирует 46800 км магистральных нефтепроводов диаметром от 420 до 1220 мм, 395 нефтеперекачивающих станций, 866 резервуаров емкостью 12,7 млн куб.м [1].

Первая страница
В данной статье описывается расчет напряженно-деформированного состояния линейно упругой пластины с отверстием. Рассматривается стационарное состояние квадратной пластины со стороной 4 м и отверстием посередине радиусом 0,5 м. К левой и правой сторонам приложены растягивающие напряжения σ=10 кПа, поскольку задача имеет две оси симметрии, то рассматривается только четвертая часть. Схема показана на рис. 3.1. Расчетная область заштрихована.

Рисунок 3.1. Геометрия плоскости с отверстием.

В данном задаче предполагаем что напряжения действующие в направлении z являются несущественными. Следовательно задача становится двумерной. Постановка задачи в плоском напряженном состоянии применяется при исследовании тонких пластин. Есть точное решение для бесконечной полуплоскости с отверстием. Выражение для напряжений по оси х имеет вид:

Рисунок 3.1. Геометрия плоскости с отверстием.
Результаты моделирования будут сравниваться с этим решением.

Сетка

Расчетная сетка состоит из пяти блоков. Как уже отмечалось в задаче о каверне. Сетка всегда является 3хмерной. Разметка приведена на рис. 3.2:

Пользователю необходимо перейти в папку plateHole в каталоге $FOAM_RUN/tutorials/stressAnalysis/solidDisplacementFoam и открыть файл constant/polyMesh/blockMeshDict. Он содержит следующий текст:
17 convertToMeters 1;
18
19 vertices
20 (
21 (0.5 0 0)
22 (1 0 0)
23 (2 0 0)
24 (2 0.707107 0)
25 (0.707107 0.707107 0)
26 (0.353553 0.353553 0)
27 (2 2 0)
28 (0.707107 2 0)
29 (0 2 0)
30 (0 1 0)
31 (0 0.5 0)
32 (0.5 0 0.5)
33 (1 0 0.5)
34 (2 0 0.5)
35 (2 0.707107 0.5)
36 (0.707107 0.707107 0.5)
37 (0.353553 0.353553 0.5)
38 (2 2 0.5)
39 (0.707107 2 0.5)
40 (0 2 0.5)
41 (0 1 0.5)
42 (0 0.5 0.5)
43 );

Пять расчетных блоков.
45 blocks
46 (
47 hex (5 4 9 10 16 15 20 21) (10 10 1) simpleGrading (1 1 1)
48 hex (0 1 4 5 11 12 15 16) (10 10 1) simpleGrading (1 1 1)
49 hex (1 2 3 4 12 13 14 15) (20 10 1) simpleGrading (1 1 1)
50 hex (4 3 6 7 15 14 17 18) (20 20 1) simpleGrading (1 1 1)
51 hex (9 4 7 8 20 15 18 19) (10 20 1) simpleGrading (1 1 1)
52 );
53
54 edges
55 (
56 arc 0 5 (0.469846 0.17101 0)
57 arc 5 10 (0.17101 0.469846 0)
58 arc 1 4 (0.939693 0.34202 0)
59 arc 4 9 (0.34202 0.939693 0)
60 arc 11 16 (0.469846 0.17101 0.5)
61 arc 16 21 (0.17101 0.469846 0.5)
62 arc 12 15 (0.939693 0.34202 0.5)
63 arc 15 20 (0.34202 0.939693 0.5)
64 );
65
66 patches
67 (
68 symmetryPlane left
69 (
70 (8 9 20 19)
71 (9 10 21 20)
72 )
73 patch right
74 (
75 (2 3 14 13)
76 (3 6 17 14)
77 )
78 symmetryPlane down
79 (
80 (0 1 12 11)
81 (1 2 13 12)
82 )
83 patch up
84 (
85 (7 8 19 18)
86 (6 7 18 17)
87 )
88 patch hole
89 (
90 (10 5 16 21)
91 (5 0 11 16)
92 )
93 empty frontAndBack
94 (
95 (10 9 4 5)
96 (5 4 1 0)
97 (1 4 3 2)
98 (4 7 6 3)
99 (4 9 8 7)
100 (21 16 15 20)
101 (16 11 12 15)
102 (12 13 14 15)
103 (15 14 17 18)
104 (15 18 19 20)
105 )
106 );
107
108 mergePatchPairs
109 (
110 );
В данной задаче мы должны использовать ячейки с изогнутыми краями. Они определяются после ключевого слова edges. Каждый элемент начинается с типа кривой, описанной позднее. Далее записываются координаты вершин дуги и координаты точек, через которые она проходит. В нашем случае это координаты середины дуги.
Не все расчетные блоки имеют одинаковую ориентацию. Как видно из рис. 3.2. направление для блока 0 отлично от ориентации для блока 4. То есть при задании сетки необходимо чтобы узлы ячеек на краях блоков совпадали.
Определяются 6 частей (патчей). По одному на каждую сторону плоскости. Одну на отверстие и одну на переднюю и заднюю сторону плоскости. Левая и нижняя части имеют симметричные плоскости, определенные при помощи команды symmetryPlane. Часть frontAndBack представляет собой плоскость, которая игнорируется в 2мерном случае. Типы расчетных частей и их геометрия подробно будут описаны позднее. Командой blockMesh формируется сетка.

Граничные и начальные условия

При отсутствии температурных напряжений начальные условия задаются в файле D.

17 dimensions [0 1 0 0 0 0 0];
18
19 internalField uniform (0 0 0);
20
21 boundaryField
22 {
23 left
24 {
25 type symmetryPlane;
26 }
27 right
28 {
29 type tractionDisplacement;
30 traction uniform ( 10000 0 0 );
31 pressure uniform 0;
32 value uniform (0 0 0);
33 }
34 down
35 {
36 type symmetryPlane;
37 }
38 up
39 {
40 type tractionDisplacement;
41 traction uniform ( 0 0 0 );
42 pressure uniform 0;
43 value uniform (0 0 0);
44 }
45 hole
46 {
47 type tractionDisplacement;
48 traction uniform ( 0 0 0 );
49 pressure uniform 0;
50 value uniform (0 0 0);
51 }
52 frontAndBack
53 {
54 type empty;
55 }
56 }

На левую и нижнюю расчетные части ставятся симметричные граничные условия при помощи команды symmetryPlane. Аналогично патч frontAndBack объявляется пустым. Остальные граничные условия описываются типом tractionDisplacement. Этот тип означает, что давление, нормальное к поверхности раздела, определяется как обычное напряжение взятое с отрицательным знаком. Поскольку напряжения на границе равны нулю, то и давление равно нулю. Для правого патча напряжение задаем следующим образом: (1е4 0 0) давление здесь также равно 0.

Физические свойства

Физические свойства устанавливаются в файле mechanicalProperties в папке constant. Физические константы для данной задачи описаны в таблице 3.1.
Таблица 3.1.

Тепловые свойства

Температурное поле T представлено в файле T в папке 0. Чтобы решить задачу с использованием температурных напряжений, необходимо установить переключатель thermalStress и положение on, в файле thermalProperties в папке constant. В таблице 3.2. показаны тепловые параметры стали.
Таблица 3.2. Тепловые характеристики стали

По умолчанию тепловое влияние не учитывается в расчете.
Управление выводом информации

Параметры выходной информации (шаг по времени, начальное/ конечное время и тд) описано в файле controlDict. Шаг по времени устанавливаем равным 1. Интервал записи в файл равным 20.

Дискретизация

Поскольку процесс стационарный, то пользователь должен указать команду SteadyState в файле fvSchemes. Уравнение движения, записанное в линейно-упругих напряжениях, включает в себя несколько членов, содержащих перемещения. Как правило, в методе конечных объемов используется метод Гаусса, которых во многих случаях дает хорошее приближение. В данном случае используется метод наименьших квадратов. Поэтому в файле fvSchemes пользователю необходимо в строке drad(T) указать командное слово leastSquares. Таким образом должна получиться следующая запись файла fvSchemes:

18 d2dt2Schemes
19 {
20 default steadyState;
21 }
22
23 gradSchemes
24 {
25 default leastSquares;
26 grad(D) leastSquares;
27 grad(T) leastSquares;
28 }
29
30 divSchemes
31 {
32 default none;
33 div(sigmaD) Gauss linear;
34 }
35
36 laplacianSchemes
37 {
38 default none;
39 laplacian(DD,D) Gauss linear corrected;
40 laplacian(DT,T) Gauss linear corrected;
41 }
42
43 interpolationSchemes
44 {
45 default linear;
46 }
47
48 snGradSchemes
49 {
50 default none;
51 }
52
53 fluxRequired
54 {
55 default no;
56 D yes;
57 T no;
58 }

Если открыть файл fvSolution, то видно, что решателем для полей D и T является GAMG. Здесь также появляются команды tolerance и relTol. Пока не очень понятно, что они означают. Кроме того появляется пункт:

45 stressAnalysis
46 {
47 compactNormalStress yes;
48 nCorrectors 1;
49 D 1e-06;
50 }

Он содержит команду nCorrectors. Этой командой определяется число циклов для расчета системы уравнений. Поскольку шаг по времени равен 1, то значение в данном случае также равно единице. Число D, определяется сходимость решения. То есть при достижении данного значения расчет для одной итерации прекращается.

Запуск расчета

Запуск расчета производится командами:
cd $FOAM_RUN/tutorials/stressAnalysis/solidDisplacementFoam/plateHole
solidDisplacementFoam > log &

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

Просмотр результата

Решатель solidDisplacementFoam выводит результаты в виде симметричного тензорного поля sigma. Просмотреть результаты можно утилитой paraFoam. Прмиер расчета показан на рис.3.3.

Рисунок 3.3. Пример расчета НДС пластины с отверстием

Сравним результат с точным решением. Необходимые данные можно получить при помощи утилиты sampleDict в папке system. Необходимые данные находятся на промежутке от (0.0, 0.5, 0.25) до (0.0, 2.0, 0.25). Сравнение графиков можно провести в GnuPlot при помощи команды:

plot [0.5:2] [0:] ’sets/100/leftPatch_sigmaxx.xy’,
1e4*(1+(0.125/(x**2))+(0.09375/(x**4)))

«sets/100/leftPatch_sigmaxx.xy» содержит расчетные данные. Пример приведен ниже:

Рисунок 3.4. Сравнение расчета и точного решения.

Упражнения

1. Увеличить число узлов сетки в 2 раза
2. Используя неравномерное разбиение по граням расчетных частей увеличить число узлов вблизи отверстия в 2 раза. Разбиение провести таким образом, чтобы отношения между соседними ячейками не превышало 1,1.
3. Оцените ошибку приближения к точному решению при изменении разбиения и размеров области (при неизменном размере отверстия).

Решения задачи в виде исполнимых файлов и описаний их приведены в моем блоге.

По всем вопросам обращаться на e-mail: serg@imech.anrb.ru

Первая страница

В данной статье описывается расчет течения в каверне. Описание данного расчета приводится в главе 2 руководства пользователя OpenFOAM 1.7.1. В каждой статье будет описываться задача: физическая постановка задачи, геометрия расчетной области, начальные и граничные условия, описание решаемых уравнений.
Все примеры находятся в обучающей директории OpenFOAM, которую можно создать на основе инсталляции программы следующими командами:
mkdir –p $FOAM_RUN
где $FOAM_RUN – наименование папки, где будут храниться обучающие программы. Например, это может быть папка LinuxOpenFOAM-1.7.1 в директории OpenFOAM-1.7.1
Обучающие программы можно затем скопировать в эту папку с помощью команды:
cp –r $FOAM_TUTORIALS/* $FOAM_RUN

Течение несжимаемой жидкости в каверне.

Геометрия расчетной области показана на рис.2.1

Рисунок 2.1: Геометрия расчетной области задачи
Область состоит из квадрата со стороной 0.1 м, 3 границы которого неподвижны, верхняя граница движется со скоростью Ux=1 м/с. Изучим возможные физические и расчетные случаи решения данной задачи.
Пользователь должен перейти в каталог cavity, командой:
cd $FOAM _RUN/tutorials/incompressible/icoFoam/cavity
Исходные данные для каждой задачи в OpenFOAM хранятся в различных файлах. Для редактирования этих файлов можно использовать любой текстовый редактор, например Emacs, VI, Gedit, Кейт, Nedit и т.д..
В OpenFOAM используется 3мерная декартова система координат. Чтобы решать задачу в 2х измерениях необходимо задать «пустые» граничные условия в направлении z. Будем решать задачу в дверной постановке. Область делится равномерной сеткой размерности 20х20 Структура разбиения показана на рисунке 2.2:

Рисунок 2.2: Структура сетки.

Файл содержащий узлы сетки генерируется командой blockMesh. Входные параметры для команды находятся в папке Polymesh в файле blockMeshDict. Которые состоит из следующих разделов:
Заголовок файла:

1 /*--------------------------------*- C++ -*----------------------------------*\
2 | ========= | |
3 | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
4 | \\ / O peration | Version: 1.7.1 |
5 | \\ / A nd | Web: www.OpenFOAM.org |
6 | \\/ M anipulation | |
7 \*---------------------------------------------------------------------------*/
8 FoamFile
9 {
10 version 2.0;
11 format ascii;
12 class dictionary;
13 object blockMeshDict;
14 }
15 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

Коэффициент масштабирования
17 convertToMeters 0.1;

Список координат вершин исследуемой области. На рис.2.2. они обозначены числами 0-7.
19 vertices
20 (
21 (0 0 0) // Вершина с номером 0
22 (1 0 0) // Вершина с номером 1
23 (1 1 0) // Вершина с номером 2
24 (0 1 0) // Вершина с номером 3
25 (0 0 0.1) // Вершина с номером 4
26 (1 0 0.1) // Вершина с номером 5
27 (1 1 0.1) // Вершина с номером 6
28 (0 1 0.1) // Вершина с номером 7
29 );

Список блоков, из которых состоит расчетная область. В данной задаче рассматриваем только один блок. Запись состоит из 3 частей. Первая часть состоит из номеров вершин данного блока. Эта часть всегда является шестигранником. Порядок этих номеров показан на рис. 2.2. Вторая запись показывает сколько точек разбиения задается в каждом направлении (направления x, y и z). В частности данная запись задает размерность задачи. Для более точного решения данной задачи можно заменить эту запись на следующую: (40 40 1). При этом шаг по времени должен быть уменьшен в 2 раза (см.ниже). Третья запись задает коэффициент расширения ячейки для каждого направления в блоке. Есть два вида расширения: упрощенный simpleGrading с 3мя параметрами (1, 2, 3), и сложный edgeGrading с параметрами расширения для каждого ребра блока. Значение определяет отношение между начальной и конечной ячейкой на грани. То есть отношение первого интервал к последнему будет равно значению описанному в simpleGrading или edgeGrading. Например, записи simpleGrading (1 2 3) и edgeGrading (1 1 1 1 2 2 2 2 3 3 3 3) будут равнозначными.

31 blocks
32 (
33 hex (0 1 2 3 4 5 6 7) (20 20 1) simpleGrading (1 1 1)
34 );
Края блоков. По умолчанию каждый край считается прямым. Поэтому в данном случае эта область пустая.
36 edges
37 (
38 );

В секции patches приведены участки сетки на которых задаются особые условия (граничные, начальные и тд.). В названии используется зарезервированное слово (patch, symmetryPlane, empty, wedge, cyclic, wall, processor) название участка сетки (патча) и его координаты. Описание зарезервированных слов будет приведено позже. Слово wall используется для описания стены. В данном случае movingWall содержит координаты движущейся стенки, fixedWalls координаты 3х неподвижных стенок. Координаты должны быть записаны в таком порядке, чтобы при направлении взгляда изнутри блока порядок вершин должен быть записан по часовой стрелке. Первая координата – любая.

40 patches
41 (
42 wall movingWall
43 (
44 (3 7 6 2)
45 )
46 wall fixedWalls
47 (
48 (0 4 7 3)
49 (2 6 5 1)
50 (1 5 4 0)
51 )

Здесь задаем пустые блоки. Слово empty означает, что на таком блоке не должно быть никаких условий. В данном случае – это передняя и задняя стенки. Порядок записи вершин такой же как и в предыдущем случае.

52 empty frontAndBack
53 (
54 (0 3 2 1)
55 (4 5 6 7)
56 )
57 );

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

59 mergePatchPairs
60 (
61 );
62
63 // ************************************************************************* //

Граничные и начальные условия

По умолчанию программой OpenFOAM считается, что начальные условия начинаются с времени t=0. Поэтому файлы с начальными условиями хранятся в папке «0». Для нашей задачи папка 0 содержит два файла: давление «p» и скорость «U». Рассмотрим файл «p»:
Определяет размерность переменной. Размерность давления - м2/с2. Полное описание данного пункта будет приведено ниже.

17 dimensions [0 2 -2 0 0 0 0];

Задание поля давления. uniform - поле давления однородно и равно 0 (в данном случае). В данном случае рассматривается изменение давления. nonuniform – значения задаются в каждой точке. Полное описание данного пункта будет приведено ниже.

19 internalField uniform 0;

Здесь описываются все граничные условия для переменных приведённых в файле blockMeshDict. Данный раздел также имеет сложный характер и будет описан ниже. Для рассматриваемой задачи граничные условия описываем равенством нулю градиента давления – zeroGradient.
21 boundaryField
22 {
23 movingWall
24 {
25 type zeroGradient;
26 }
27
28 fixedWalls
29 {
30 type zeroGradient;
31 }

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

33 frontAndBack
34 {
35 type empty;
36 }
37 }

Аналогично описываются значения для скорости в файле «U»:

17 dimensions [0 1 -1 0 0 0 0];

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

19 internalField uniform (0 0 0);

Задание граничных условий отличается тем, что мы задаем значение равное 1 м/с по переменной х. Правильное описание приведено в секции movingWall. На всех остальных стенках начальная скорость равна 0.
21 boundaryField
22 {
23 movingWall
24 {
25 type fixedValue;
26 value uniform (1 0 0);
27 }
28
29 fixedWalls
30 {
31 type fixedValue;
32 value uniform (0 0 0);
33 }
34
35 frontAndBack
36 {
37 type empty;
38 }
39 }

Физические свойства.

Все физические свойства переменных приведены в начале файлов. В файле transportProperties для рассматриваемой задачи приведено значение для кинематической вязкости и её размерность:

18 nu nu [ 0 2 -1 0 0 0 0 ] 0.01;

Кинематическая вязкость записывается фонетическим звучанием греческой буквы «nu». Для течения несжимаемой жидкости в каверне используется следующая формула:

где d и U характерная длина стенки каверны и скорость. В данном случае d=0,1 м, |U|=1 м/с.
Данные, связанные с контролем шага по пространству и времени считываются из файла controlDict. Он находится в папке system. В примере используется начальное время равное 0, значит начальные данные OpenFOAM будет считывать из папки под наименованием 0. Начальное время устанавливается в секции startTime. Для окончания счета необходимо в секции stopAt записать слово endTime. В данной задаче endTime равно 0,5. Теперь необходимо установить шаг по времени deltaT. Для устойчивого счета необходимо выполнение следующего условия:

где dt - шаг по времени, dx - размер ячейки. Выбираем dt из условия, что максимум Co соответствует максимальной скорости и минимальной ячейке. Из этих условий следует, что dt =0,005 с. Стоит заметить, что в данной задаче выбор промежутка времени определяется пользователем. Для записи поля скоростей и давления в определенные промежутки времени в секцию writeControl записываем слово timeStep, которое означает, что запись будет происходить через равное число итераций, описанное в секции writeInterval.

Дискретизация и настройка решения

Пользователь указывает выбор схемы дискретизации в файле fvSchemes в папке system.
Подробнее данный раздел будет описан позже.

Просмотр сетки

В некоторых случаях стоит перед расчетом просмотреть исследуемую область. Разработчики OpenFOAM рекомендуют для этих целей использовать свободное программное обеспечение ParaView. Данная программа запускается из директории cavity (для текущей задачи) командой ParaFoam. Работа с paraView будет описана позже.

Запуск расчета

Как любое UNIX приложение OpenFOAM может быть запущено двумя способами. Либо в каталоге cavity командой:

icoFoam

либо в командной строке с использованием аргумента -case:

icoFoam -case $FOAM_RUN/tutorials/incompressible/icoFoam/cavity

Параметры расчета высвечиваются в окне терминала. Чтобы записать процесс расчета необходимо выполнить расчет при помощи следующей команды (из директории cavity)

icoFoam > log &
cat log

Обработка полученного решения

Расчетная задача в ПС OpenFOAM представляется как упорядоченный набор файлов и каталогов, размещенных в директории файловой системы под некоторым именем. Корневая структура расчетной задачи (относительно каталога задачи) имеет вид, показанный на рис. 2.3..

Рисунок 2.3.Структура расчетных файлов OpenFOAM
После расчета результаты можно просмотреть командой paraFoam в окне программы paraView.

Создание сложной сетки

Если в некоторой области значения скорости или давления сильно отличаются в ближайших точках то необходимо использовать более сложные виды расчетных сеток. В нашем случае в верхних крайних точках и в области около них значения давления жидкости гораздо больше чем в центре каверны. Зачастую возникают такие задачи, где в некоторых областях достаточно знать лишь очень грубое приближение. Например, задача о распространении трещины гидроразрыва. В окрестности вершины трещины необходимо знать очень точные значения напряжений, на расстояниях много больше диаметра трещины можно использовать более грубую сетку.
Исследуем задачу течения жидкости в каверне для сложной сетки. Поскольку изменения в исходных файлах довольно существенны, пример решения такой задачи находится в папке CavityGrade в стандарной поставке OpenFOAM.
Создадим сетку как показано на рис. 2.4:

Рисунок.2.4. Структура сетки для задачи течения несжимаемой жидкости в каверне.

Файл blockMeshDict будет иметь следующий вид:

17 convertToMeters 0.1;
18
19 vertices
20 (
21 (0 0 0)
22 (0.5 0 0)
23 (1 0 0)
24 (0 0.5 0)
25 (0.5 0.5 0)
26 (1 0.5 0)
27 (0 1 0)
28 (0.5 1 0)
29 (1 1 0)
30 (0 0 0.1)
31 (0.5 0 0.1)
32 (1 0 0.1)
33 (0 0.5 0.1)
34 (0.5 0.5 0.1)
35 (1 0.5 0.1)
36 (0 1 0.1)
37 (0.5 1 0.1)
38 (1 1 0.1)
39 );

Здесь уже 4 блока. Ребро каждого разбивается на 10 отрезков длины которых меняются с коэффициентами, описанными командой simpleGrading.

41 blocks
42 (
43 hex (0 1 4 3 9 10 13 12) (10 10 1) simpleGrading (2 2 1)
44 hex (1 2 5 4 10 11 14 13) (10 10 1) simpleGrading (0.5 2 1)
45 hex (3 4 7 6 12 13 16 15) (10 10 1) simpleGrading (2 0.5 1)
46 hex (4 5 8 7 13 14 17 16) (10 10 1) simpleGrading (0.5 0.5 1)
47 );
48
49 edges
50 (
51 );
52
53 patches
54 (

Поскольку верхняя стенка теперь состоит из 2 блоков, то секция movingWall будет состоять из 2х записей,

55 wall movingWall
56 (
57 (6 15 16 7)
58 (7 16 17 8)
59 )

то же относится и к секции fixedWalls.

60 wall fixedWalls
61 (
62 (3 12 15 6)
63 (0 9 12 3)
64 (0 1 10 9)
65 (1 2 11 10)
66 (2 5 14 11)
67 (5 8 17 14)
68 )

В 4 раза увеличивается количество пустых граничных условий.

69 empty frontAndBack
70 (
71 (0 3 4 1)
72 (1 4 5 2)
73 (3 6 7 4)
74 (4 7 8 5)
75 (9 10 13 12)
76 (10 11 14 13)
77 (12 13 16 15)
78 (13 14 17 16)
79 )
80 );
81
82 mergePatchPairs
83 (
84 );

Вычисление шага по времени

Максимальная скорость движения жидкости и минимальные ячейки находятся около верхней грани каверны. Использование команды simpleGrading позволяет задавать ячейки с использованием геометрической прогрессии. Если длина ребра блока равна l, а число точек на ребре n, с отношением между первым отрезком и последним R, то размер самой меленькой ячейки определяется по следующей формуле:

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

В рассматриваемом случае dxs=0,00345 м. Следовательно, временной интервал равен dt=0,0025 с. Поэтому значение writeInterval (файл controlDict в папке cavityGrade/system/) устанавливаем равным 40, чтобы расчетные значения записывались через каждые 0,1 с. Теперь необходимо установить начальное и конечное время работы программы и произвести расчет командами:
blockMesh
icoFoam
после расчета можно просмотреть результат командой:
paraFoam
Стоит заметить, что в OpenFOAM совпадение узлов расчетной сетки на границах блоков полностью определяется пользователем. Если узлы не будут совпадать, то расчет будет неверным.

Увеличение числа Рейнольса.

При увеличении числа Рейнольдса до 100, схождение к стационарному решению будет происходить гораздо медленнее. Для решения данной задачи необходимо создать папку cavityHighRe. В стандартной поставке OpenFOAM в папке icoFoam находится исполнимый файл Allrun. После работы этого файла в папке появляются все директории описанные в данном руководстве. Папка cavityHighRe создается при помощи следующих команд:

cd $FOAM_RUN/tutorials/incompressible/icoFoam
cp -r cavity cavityHighRe

Войдите в файл transportProperties в папке constant. Поскольку мы увеличиваем число Рейнольдса в 10 раз, то необходимо уменьшить вязкость в 10 раз. Так как стационарных режим возникает позже, то необходимо расчетное время увеличить в 2 раза (в файле cavityHighRe/system/controlDict в секции endTime необходимо поставить число 2).
Теперь запустите расчет командой icoFoam. При выполнении работы в фоновом режиме, могут быть полезны следующие команды UNIX:

nohup – программа остается в работе, даже если пользователь закрыл терминал.
nice - изменяет приоритет процесса. -20 – наивысший приоритет, 19 – низший.

Это полезно, например, если пользователь хочет провести расчет на удаленной машине и установить ему низкий приоритет. Пример выполнения команд:

cd $FOAM_RUN/tutorials/incompressible/icoFoam/cavityHighRe
nohup nice -n 19 icoFoam > log &
cat log

Часто время установления стационарного режима гораздо меньше чем заданное пользователем. Для того чтобы программа не работала впустую в файле fvSolutions можно задать разницу между ближайшими итерациями скорости U. Обычно это значение равно 10^-6. Как только разница стала меньше чем заданное число расчет останавливается.
Разумеется бесконечно увеличивать число Рейнольса нельзя. Поскольку возникают турбулентные режимы, которые описываются другими законами. Исследование турбулентности и более сложные режимы течения можно рассчитать перейдя в папку pisoFoam.

Сложные режимы течения

Перейдите в папку $FOAM_RUN/tutorials/incompressible/pisoFoam/ras. Сгенерируйте сетку командой blockMesh. Начиная с версии 1.6. в OpenFOAM применяются специальные модели, позволяющие задать разные граничные условия на различных стенках. Некоторые из них приведены при выборе турбулентной вязкости в файле 0/nut:
18 dimensions [0 2 -1 0 0 0 0];
19
20 internalField uniform 0;
21
22 boundaryField
23 {
24 movingWall
25 {
26 type nutWallFunction;
27 value uniform 0;
28 }
29 fixedWalls
30 {
31 type nutWallFunction;
32 value uniform 0;
33 }
34 frontAndBack
35 {
36 type empty;
37 }
38 }
Здесь используется граничное условие nutWallFunction. В файлах 0/k и 0/epsilon используются граничные условия epsilonWallFunction и kqRwallFunction. Начальные значения для k и ε выбираются на основе следующих формул:


где - константа равная 0,09 для k-ε модели. В декартовых координатах k определяется по формуле скалярного произведения:

где - координаты скорости. Предположим что начальная турбулентность изотропна, т. е. и равна 5% от скорости верхней грани каверны. Тогда получаем:



До версии 1.6 для каждого крупного вихря создавался свой решатель, это приводило к дублированию кода. Теперь метод моделирования выбирается пользователем при помощи ключевого слова simulationType в файле turbulenceProperties. В папке можно открыть этот файл:

18 simulationType RASModel;

Можно использовать команды laminar, RASModel и LESModel. В данном случае выбирается модель RAS определенная в файле RASProperties в папке constant. Виды RAS моделей будут описаны позднее. Здесь kEpsilon – это стандартная k-ε модель турбулентности. В данном случае пользователь должен установить флаг on для расчета турбулентности:

turbulence on;

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

Упражнение

1. Установите входные параметры таким образом, чтобы расчет происходил быстрее, при этом не теряя точности.

Изменение геометрии расчетной области

Создадим папку cavityClipped в директории icoFoam. Для примера удалим квадрат со стороной 0,04 м из расчетной области (внизу справа). Файл BlockMeshDict будет выглядеть следующим образом:

17 convertToMeters 0.1;
18
19 vertices
20 (
21 (0 0 0)
22 (0.6 0 0)
23 (0 0.4 0)
24 (0.6 0.4 0)
25 (1 0.4 0)
26 (0 1 0)
27 (0.6 1 0)
28 (1 1 0)
29
30 (0 0 0.1)
31 (0.6 0 0.1)
32 (0 0.4 0.1)
33 (0.6 0.4 0.1)
34 (1 0.4 0.1)
35 (0 1 0.1)
36 (0.6 1 0.1)
37 (1 1 0.1)
38
39 );
40

Расчетная область состоит из 3 блоков.
41 blocks
42 (
43 hex (0 1 3 2 8 9 11 10) (12 8 1) simpleGrading (1 1 1)
44 hex (2 3 6 5 10 11 14 13) (12 12 1) simpleGrading (1 1 1)
45 hex (3 4 7 6 11 12 15 14) (8 12 1) simpleGrading (1 1 1)
46 );
47
48 edges
49 (
50 );
51
52 patches
53 (
54 wall lid
55 (
56 (5 13 14 6)
57 (6 14 15 7)
58 )
59 wall fixedWalls
60 (
61 (0 8 10 2)
62 (2 10 13 5)
63 (7 15 12 4)
64 (4 12 11 3)
65 (3 11 9 1)
66 (1 9 8 0)
67 )
68 empty frontAndBack
69 (
70 (0 2 3 1)
71 (2 5 6 3)
72 (3 6 7 4)
73 (8 9 11 10)
74 (10 11 14 13)
75 (11 12 15 14)
76 )
77 );
78
79 mergePatchPairs
80 (
81 );

Поскольку в поставке OpenFOAM учебные файлы в данном случае проводят расчет со времени 0,5 c. То пользователю необходимо записать начальные данные в папку 0,5 предварительно создав её. Сделать это можно следующими командами:

cd $FOAM_RUN/tutorials/incompressible/icoFoam/cavityClipped
cp -r 0 0.5

После всех изменений можно провести расчет командой icoFoam.

Упражнение

1. Проделайте то же самое, только удалите квадрат сверху слева.
2. Проведите расчет для более сложной области, в виде квадрата удалённого в середине расчетной области.
3. Проведите расчет для полостей в виде прямоугольников со сторонами 0,04 и 0,08 м удаленными слева и справа внизу расчетной области.

Решения упражнений в виде исполнимых файлов и описание их будут с течением времени приведены в моем блоге.

По всем вопросам обращаться на e-mail: serg@imech.anrb.ru

Первая страница
Рассмотрим следующие задачи:

1. Установка OpenFOAM 1.7.1. через Интернет в UBUNTU
2. Инсталляция из исходных файлов
3. Инсталляция в дистрибутивах, отличных от UBUNTU

Установка OpenFOAM 1.7.1. через Интернет в UBUNTU

Для корректной инсталляции последней версии OpenFOAM необходимы следующие версии UBUNTU:

• 10.04 LTS, codename lucid
• 10.10, codename maverick

В терминале заходим под паролем администратора (root) при помощи команды «sudo -i». Набираем следующую команду:

sudo sh -c "echo deb http://www.openfoam.com/download/ubuntu main >> /etc/apt/sources.list"

здесь вместо надо набрать либо слово lucid, либо слово maverick, в зависимости от версии Ubuntu. Например в случае UBUNTU 10.04 LTS, codename lucid данная команда будет выглядеть так:

sudo sh -c "echo deb http://www.openfoam.com/download/ubuntu lucid main >> /etc/apt/sources.list"

Данную команду необходимо выполнить только один раз.
Обновляем установочные пакеты:
sudo apt-get update
Инсталлируем новую версию. На текущий момент – это версия 1.7.1.:
sudo apt-get install openfoam171
Инсталлируем Paraview (на текущий момент – версия 3.8.1):
sudo apt-get install paraviewopenfoam381
Теперь OpenFOAM-1.7.1 и Paraview-3.8.1 находятся в директории /opt.
Для полного завершения установки необходимо настроить OpenFOAM при помощи следующих команд:
открыть файл .bashrc в домашней директории пользователя. Либо набрать в терминале:
gedit ~/.bashrc
В самом низу файла добавить следующую строку:
. /opt/openfoam171/etc/bashrc
Если там уже есть подобная строка, то её необходимо закомментировать символом #.
После этого сохраняем документ и выполняем команду
. ~/.bashrc
OpenFOAM установлен и готов к работе.

Инсталляция из исходных файлов

Необходимо скачать следующие файлы:
OpenFOAM - OpenFOAM-1.7.1.gtgz
Third-Party - ThirdParty-1.7.1.gtgz

Затем необходимо выбрать папку, куда будут распакованы данные файлы, эта папка станет инсталляционной для OpenFOAM. Если пользователем программы будет только один человек, либо пользователь не хочет устанавливать продукт под правами администратора рекомендуется установить систему в домашнюю папку. Если наоборот, то можно установить в одну из следующих папок: /usr/local/OpenFOAM, /opt/OpenFOAM, либо просто /opt.
После создания папки скопируйте туда 2 сохраненных файла и распакуйте их следующими командами:

• tar xzf OpenFOAM-1.7.1.gtgz
• tar xzf ThirdParty-1.7.1.gtgz

Файлы распакуются в папки OpenFOAM-1.7.1 и ThirdParty-1.7.1.

Системные требования.

OpenFOAM-1.7.1 и ThirdParty-1.7.1 тестировались на Linux дистрибутивах Ubuntu 10.04 и OpenSuSE 11.1, 11.2 и 11.3. ParaView-3.8.0 требует наличие cmake-2.6.4 или выше и Qt-4.5 или выше. Чтобы протестировать свою систему необходимо выполнить скрипт foamSystemCheck, описанный в разделе “Тестирование системы ” (ниже).

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

Файлы переменных среду окружения содержатся в папке OpenFOAM-1.7.1/etc. Например, если система установлена в домашнюю папку $HOME/OpenFOAM, то путь к переменным будет следующий:
• $HOME/OpenFOAM/OpenFOAM-1.7.1/etc
либо
при работе с bash или ksh, необходимо добавить следующую строку в конец файла .bashrc:
• . $HOME/OpenFOAM/OpenFOAM-1.7.1/etc/bashrc
затем обновить командой “. $HOME/.bashrc” в текущем окне терминала
либо
при работе с tcsh или csh, в папке etc/cshrc есть файл .cshrc в конец его надо добавить:
• source $HOME/OpenFOAM/OpenFOAM-1.7.1/etc/cshrc
затем обновить командой “source $HOME/.cshrc” в текущем окне терминала.

Тестирование системы

Для проверки пригодности системы необходимо выполнить скрипт foamSystemCheck (он находится в каталоге OpenFOAM-1.7.1/bin). Если какой-либо нужное обеспечение отсутствует, или нуждается в обновлении до новой версии, установите необходимое программное обеспечение, прежде чем приступить к инсталляции.
В частности: на вашем компьютере должна быть последняя версия GCC (достаточно версии-4.3., рекомендуется версия 4.4.). Последняя версия доступна на сайте http://gcc.gnu.org/ .
Чтобы инсталлировать OpenFOAM вам также понадобится программное обеспечение Flex.
Инсталляция из исходников
Перейти в самую верхнюю папку $WM_PROJECT_DIR и выполнить скрипт ./Allwmake. Если возникнут трудности, свяжитесь с поддержкой .

Инсталляция в дистрибутивах, отличных от UBUNTU

При установке на другие дистрибутивы, можно руководствоваться следующими источниками:

1. Установка OpenFOAM 1.5 в Ubuntu 8.10
2. Как установить OpenFOAM на Ubuntu 9.04?
3. Установка OpenFOAM на ALTLinux 4 Desktop
4. Установка OpenFOAM 1.7.0 на OpenSUSE 11.3 (с нуля) (PDF)
5. OpenFOAM 1.7.x on Windows 64-bit with Native MPI
6. Дистрибутив CAELinux

e-mail: serg@imech.anrb.ru

Первая страница

Введение

На текущий момент существует огромное количество программ математического моделирования (ANSYS, LabVIEW, MATLAB, Maple, Mathematica и тд). В основном это коммерческое программное обеспечение. Поскольку эти программы стоят достаточно дорого, для них есть фирмы, обеспечивающие поддержку и разработку документации, группа клиентов. Кроме них существует также программное обеспечение с открытым исходным кодом (OpenFOAM, SALOME, KCalc, Scilab, R и тд). В основном, это программы, разработанные под операционную систему Linux. Одним из минусов открытого программного обеспечения является разнородность групп поддержки, практически полное отсутствие документации, в том числе и русской. Данное издание предназначено для обучения работе в OpenFOAM. Предполагается, что читатель знаком с основами дифференциального исчисления в частных производных, основами объектно-ориентированного программирования на языке C++, владеет основами программирования для CUDA

OpenFOAM (Open Source Field Operation And Manipulation) — открытая интегрируемая платформа для численного моделирования задач механики сплошных сред. Предназначена для решения задач в области механики сплошных сред.
По данным сайта http://www.openfoam.com/ 10 декабря 2004 года код FOAM (Field Operation And Manipulation) был открыт под лицензией GNU GPL. Таким образом, данная дата считается первым релизом открытого программного продукта OpenFOAM. FOAM остался закрытым и продолжает развиваться параллельно с OpenFOAM. На текущий момент доступна для скачивания версия 1.7.1. анонсированная 26 августа 2010 года.
OpenFOAM позволяет решать следующие задачи:

• Прочностные расчеты
• Гидродинамика ньютоновских и неньютоновских вязких жидкостей как в несжимаемом, так и сжимаемом приближении с учётом конвективного теплообмена и действием сил гравитации. Для моделирования турбулентных течений возможно использование RANS-моделей, LES и DNS методов. Возможно решение дозвуковых, околозвуковых и сверхзвуковых задач;
• Задачи теплопроводности в твёрдом теле; (
• Многофазные задачи, в том числе с описанием химических реакций компонент потока;
• Задачи, связанные с деформацией расчётной сетки;
• Сопряжённые задачи;
• Некоторые другие задачи, при математической постановке которых требуется решение ДУ в ЧП в условиях сложной геометрии среды;
• Распараллеливание расчёта как в кластерных, так и многопроцессорных системах



Организация OpenFOAM

Условно код OpenFOAM можно разделить на следующие части (как и любую другую сложную программу):
1. Ядро (core) – набор файлов, определяющих основные функциональные возможности кода
2. Библиотеки (solvers) – содержат множество моделей механики сплошных сред.
3. Исполняемые файлы (utilities)– утилиты для решения задач МСС, подготовки и обработки данных

• Утилиты для работы с расчетной сеткой
• Диагностические утилиты
• Утилиты для задания начальных и граничных условий
• Решатели
• Утилиты для обработки полученных данных
• Системные утилиты

Рабочим языком OpenFOAM является язык C++. Данный язык выбран для удобства представления математических операторов и выражения в удобной форме до компиляции в исполняемый файл. Сам продукт состоит из инструментов для решения систем дифференциальных уравнений в частных производных (solvers), утилит (utilities) и библиотек (libraries). При решении задач используется конечно-элементный метод.
В качестве примера обычно приводится уравнение сохранения количества движения для ньютоновской несжимаемой жидкости без действия массовых сил:

которое может быть представлено в виде:
solve
(
fvm::ddt(rho, U)
+ fvm::div(rho, U, U)
- fvm::laplacian(mu, U)
==
- fvc::grad(p)
);

которое записывается следующим образом:
solve
(
fvm::ddt(k)
+ fvm::div(phi, k)
- fvm::laplacian(nu() + nut, k)
== nut*magSqr(symm(fvc::grad(U)))
- fvm::Sp(epsilon/k, k)
);
Другие примеры:
Уравнение неразрывности:
solve(
fvm::ddt(rho)
+ fvc::div(phi)
);
Уравнение сохранения количества движения:
fvVectorMatrix UEqn
(
fvm::ddt(rho, U)
+ fvm::div(phi, U)
- fvm::laplacian(mu, U)
);
solve(UEqn == -fvc::grad(p));
Уравнение сохранения энергии:
solve
(
fvm::ddt(rho, e)
+ fvm::div(phi, e)
- fvm::laplacian(mu, e)
==
- p*fvc::div(phi/fvc::interpolate(rho))
+ mu*magSqr(symm(fvc::grad(U)))
);
T = e/Cv;
Официальный сайт в состав поставки включает Paraview 3.8.1 - программа для графической визуализации данных.
Часто для моделирования исходной сетки используется программный продукт SALOME.
Компании и академические институты, которые используют OpenFOAM: ABB Corporate Research, Audi, Airbus, Bayer, CD-adapco group, Danone, Daimler, Hitachi, Mitsubishi, Obayashi, SKF, Shell, Toyota, Tokyo Gas, Volkswagen, Electrolux, Fluid Technology, Icon-CG, The Technology Partnership plc, Imperial College London, Kings College London, Chalmers University, University of Exeter, University of Guelph, Hirosaki Universty, Tokyo Institute of Technology.

При подготовке данной статьи использовались материалы из следующих источников:

• Официальный сайт
• Сайт обмена практическим опытом и использования свободного ПО в CAD/CAE
• http://openfoamwiki.net/index.php/Main_Page
• http://www.foamcfd.org/

e-mail: serg@imech.anrb.ru

Лабораторией МТП разработан комплекс задач по расчету оптимальных режимов работ магистральных нефтепроводов по выбранным пользователем параметрам.
Комплекс позволяет решать следующие задачи:

1. Моделирование технологических схем работы нефтепроводов.
2. Идентификация эксплуатационных параметров линейной части и насосно-силового оборудования нефтепроводов.
3. Разработка оптимальных режимов работы нефтепроводов
4. Определение объема перекачки по нефтепроводам при заданном (установленном) включенном в работу насосном оборудовании.
5. Разработка технологических карт режимов работы нефтепроводов со сбросами и подкачками.
6. Определение перевальных точек и параметров их существования и ликвидации.
7. Разработка оптимального режима работы нефтепроводов для заданного объема перекачки.

Для получения дополнительной информации о программном продукте необходимо написать письмо на email: serg@uimech.org

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

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

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

Содержание:

1. Математическая модель
2. Конечно-разностная аппроксимация системы уравнений
3. численные расчёты (сравнение с экспериментом, сравнение распространения волн в консолидированной среде и среде насыпной плотности)
4. Влияние переменной пористости на характер волнового течения
5. Влияние переменной пористости на характер волнового течения

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

Рассмотрим плоское одномерное движение двухфазной среды на основе физических законов сохранения. Нижний индекс i = 1 будем относить к параметрам жидкости, а i = 2 – к параметрам дисперсной фазы (твердых частиц). Считаем, что в процессе распространения волн давления в чистой жидкости и в пористой среде температуры фаз изменяются незначительно.

где , - приведённая плотность и массовая скорость i-й фазы, p1, – давление и пористость жидкости, , – объёмное содержание и приведённое напряжение твёрдой фазы. Напряжение интерпретируется как часть тензора напряжений твёрдой фазы , обусловленная передачей усилий через контакты между зёрнами и не зависящая от давления жидкости. F12 – межфазная сила, связанная со скоростной неравновесностью фаз.
Сила взаимодействия F12 между несущей средой и включениями обязана трём основным эффектам:

Здесь – сила трения Стокса. Возникает из-за действия вязких сил при взаимодействии между фазами. – сила, связанная с воздействием «присоединённых масс», возникающая из-за ускорения движения включения относительно несущей среды, когда в последней возникают возмущения на расстоянии порядка размера включений. И, – сила Бассэ-Буссинеска, действующая на сферу во флюиде в момент t и зависящая от всей предыстории движения.
Для указанных сил можно записать следующие соотношения:

Влияние формы частиц и структура среды учитывается в коэффициентах , и . Их значения следует выбирать при сопоставлении с реперным экспери-ментом.
Уравнение состояния вязкоупругого скелета можно представить в следующем виде:

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

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

где - начальная приведённая плотность частиц скелета, Df*, De* – скорости звука в скелете. Данные соотношения обусловлены точностью эксперимен-тального определения скорости звука в скелете по сравнению с непосредствен-ным измерением модуля упругости. С теоретической точки зрения подобное представление модулей создаёт дополнительные возможности при построении моделей, учитывающих зависимости параметров среды от пористости. Примем, что твёрдые вещества – металлы, частицы горных пород, а также жидкости при действии всестороннего сжатия подчиняются линейному закону, поэтому для материалов фаз запишем следующие уравнения состояния:

где - скорость звука в материале i – ой фазы, , - её

сжимаемость и начальная истинная плотность.
В соответствии с уравнением совместности деформаций имеет место и следующее равенство:

Таким образом, получили систему уравнений, описывающую движение одномерной двухфазной среды:

Вопросы гиперболичности, устойчивости и корректности задачи Коши применительно к системе дифференциальных уравнений двухскоростного дви-жения дисперсных сред рассматривались в работах [39-41, 73].

Конечно-разностная аппроксимация системы уравнений

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