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

Первая страница
В данной статье описывается расчет напряженно-деформированного состояния линейно упругой пластины с отверстием. Рассматривается стационарное состояние квадратной пластины со стороной 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