5 Вычисление эффективных свойств (гомогенизация) периодических структур

5.1 Аннотация и цели обучения

Предварительные знания: Основы линейной теории упругости, нотация Фойгта (Урок 2), навыки работы с глобальными матрицами жесткости, тензорная алгебра.
Освоенные численные методы: Метод асимптотической гомогенизации, концепция представительного элемента объема (RVE), генерация строго периодических сеток в Gmsh, реализация периодических граничных условий (PBC) методом множителей Лагранжа через матрицу алгебраических ограничений.
Физическая задача: Определение макроскопических эффективных упругих характеристик (тензора упругости) гетерогенной периодической сотовой структуры, применяемой в аддитивных технологиях и при создании метаматериалов.

5.2 Осовные теоретические сведения

5.2.1 Асимптотическая гомогенизация и масштабный переход

Многие современные материалы, такие как композиты, пенопласты, 3D-печатные сотовые структуры, имеют ярко выраженную внутреннюю микроструктуру. Моделировать каждую пору или каждое волокно в масштабах крупной детали (например, крыла самолета) вычислительно невозможно.

Метод асимптотической гомогенизации решает эту проблему путем разделения задачи на два масштаба:

Микромасштаб ($l$): Уровень отдельной ячейки материала, где явно моделируются все пустоты, включения и геометрические особенности. На этом уровне материалы компонент считаются сплошными с известными свойствами.
Макромасштаб ($L$): Уровень всей конструкции ($l \ll L$). На этом уровне материал считается однородным (гомогенным), но обладающим эффективными (усредненными) свойствами, вычисленными из микромасштаба.

Связующим звеном выступает представительный объемный элемент (Representative Volume Element, RVE). Для структур с периодической геометрией (например, сотовых) RVE представляет собой одну минимальную повторяющуюся ячейку.

Фундаментальным условием корректного перехода между масштабами является теорема Хилла - Манделя, которая гласит, что макроскопическая работа деформации должна быть равна среднему по объему значению микроскопической работы:

\[ \langle \mathbf{\sigma} : \mathbf{\varepsilon} \rangle = \langle \mathbf{\sigma} \rangle : \langle \mathbf{\varepsilon} \rangle \]

где скобки $\langle \cdot \rangle$ обозначают осреднение величины по объему RVE: $\langle f \rangle = \frac{1}{V}\int_V f \, dV$.

5.2.2 Осреднение и краевые задачи

На микроуровне внутри RVE (область $\Omega$) справедливы классические уравнения линейной теории упругости:

\[ \nabla \cdot \mathbf{\sigma} = 0, \quad \mathbf{\sigma} = {}^4\mathbf{C} : \mathbf{\varepsilon}(\mathbf{u}) \]

После процедуры осреднения, связь между макроскопическими (эффективными) напряжениями и деформациями принимает вид обобщенного закона Гука:

\[ \langle \mathbf{\sigma} \rangle = {}^4\mathbf{C}^{\mathrm{eff}} \langle \mathbf{\varepsilon} \rangle \]

где ${}^4\mathbf{C}^{\mathrm{eff}}$ - тензор эффективных упругих компонент четвертого порядка. Однако часто его представляют в виде матрицы используя нотацию Фойгта (см. Урок 2) и в двумерном случае матрица эффективных компонент $\mathbf{C}_v^{\mathrm{eff}}$ имеет размерность $3 \times 3$. Чтобы найти все её компоненты, нам необходимо решить три независимые краевые задачи, прикладывая к ячейке три макроскопических состояния деформации $\bar{\varepsilon}$:

Одноосное растяжение вдоль оси X: $\bar{\varepsilon}_{11} = 1$
Одноосное растяжение вдоль оси Y: $\bar{\varepsilon}_{22} = 1$
Чистый сдвиг в плоскости XY: $\bar{\varepsilon}_{12} = 1$

Прикладывая единичную макро-деформацию, осредненный вектор напряжений $\langle \mathbf{\sigma}_v \rangle$ будет в точности равен соответствующему столбцу матрицы $\mathbf{C}_v^{\mathrm{eff}}$.

5.2.3 Типы граничных условий

Существует три способа приложить макро-деформацию $\bar{\varepsilon}$ к границе ячейки:

Кинематические (KUBC): Жесткое задание перемещений на границе. Завышает жесткость.
Статические (SUBC): Задание усилий (тракций) на границе. Занижает жесткость.
Периодические (PBC): Противоположные границы ячейки могут искривляться свободно, но их перемещения строго связаны друг с другом:

\[ \mathbf{u}(\mathbf{x}^+) - \mathbf{u}(\mathbf{x}^-) = \bar{\varepsilon} \cdot (\mathbf{x}^+ - \mathbf{x}^-) \]

где $\mathbf{x}^+$ и $\mathbf{x}^-$ — координаты соответствующих точек на противоположных гранях.

Для строго периодических сред метод PBC дает точное значение эффективной жесткости. Математически доказано, что $\mathbf{C}^{SUBC} \le \mathbf{C}^{PBC} \le \mathbf{C}^{KUBC}$.

Для реализации PBC нам необходимо построить идеальную периодическую геометрию и связать степени свободы противоположных узлов. Это приведет к системе уравнений с ограничениями (алгебраическими связями). В методе конечных элементов такие задачи решаются с помощью множителей Лагранжа, где исходная матрица жесткости $\mathbf{K}$ дополняется матрицей ограничений $\mathbf{A}$:

\[ \begin{bmatrix} \mathbf{K} & \mathbf{A}^T \\ \mathbf{A} & \mathbf{0} \end{bmatrix} \begin{Bmatrix} \mathbf{u} \\ \mathbf{\lambda} \end{Bmatrix}= \begin{Bmatrix} \mathbf{0} \\ \mathbf{q} \end{Bmatrix} \]

Вектор $\mathbf{q}$ задает относительные смещения границ, определяемые макроскопической деформацией $\bar{\varepsilon}$.

5.3 Практический задача

5.3.1 Исходные данные

Смоделируем 2D-ячейку полимерного материала полученную при помощи 3D-печати (сотовая структура). Материал: $E =$ 70 ГПа, $\nu =$ 0.33.

5.3.2 Построение геометрии и периодической сетки в Gmsh

Репрезентативная ячейка сотовой структуры (Рисунок 5.1 (a)) представляет собой прямоугольник с вырезанными половинками правильных шестиугольников по краям, в зависимости от того, как определить границы репрезентативной ячейки (см. Рисунок 5.1 (b)).

Для построения геометрии ячейки, создадим вспомогательную функцию create_hexagon, которая будет создавать шестиугольник. Далее создаем прямоугольник $l_x \times l_y$ и 4 шестиугольника, которые в дальнейшем вырезаются из прямоугольника помощи булевой операции cut (см. Список 5.1).

Список 5.1: Генерация геометрии Honeycomb и наложение периодичности сетки

gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 1)
gmsh.model.add("Honeycomb_RVE")

R_pitch = 1.0       # Шаг сот 
t = 0.3             # Толщина стенки напечатанного материала

R_hole = R_pitch - t / sqrt(3)
lx = sqrt(3) * R_pitch  # Ширина RVE
ly = 3.0 * R_pitch      # Высота RVE

# Функция создания шестиугольника
function create_hexagon(x, y, R)
    pts = []
    # Поворачиваем на 30 градусов, чтобы грани были горизонтальными
    for i in 0:5
        angle = i * pi / 3 + pi / 6
        px = x + R * cos(angle)
        py = y + R * sin(angle)
        push!(pts, gmsh.model.occ.addPoint(px, py, 0))
    end
    lines = []
    for i in 1:6
        p1 = pts[i]
        p2 = pts[i == 6 ? 1 : i + 1]
        push!(lines, gmsh.model.occ.addLine(p1, p2))
    end
    cl = gmsh.model.occ.addCurveLoop(lines)
    return gmsh.model.occ.addPlaneSurface([cl])
end

rect = gmsh.model.occ.addRectangle(0, 0, 0, lx, ly)
hex_tags = [
    create_hexagon(lx / 2, 0, R_hole),
    create_hexagon(lx / 2, ly, R_hole),
    create_hexagon(0, ly / 2, R_hole),
    create_hexagon(lx, ly / 2, R_hole),
]

gmsh.model.occ.cut([(2, rect)], [(2, tag) for tag in hex_tags])
gmsh.model.occ.synchronize()

EPS = 1e-6
# Матрицы трансляции для периодической сетки
x_translation = [1 0 0 lx; 0 1 0 0; 0 0 1 0; 0 0 0 1]
y_translation = [1 0 0 0; 0 1 0 ly; 0 0 1 0; 0 0 0 1]

# Поиск парных границ
bot_left = gmsh.model.occ.get_entities_in_bounding_box(-EPS, -EPS, -EPS, lx / 2 + EPS, +EPS, +EPS, 1)
top_left = gmsh.model.occ.get_entities_in_bounding_box(-EPS, ly - EPS, -EPS, lx / 2 + EPS, ly + EPS, +EPS, 1)
bot_right = gmsh.model.occ.get_entities_in_bounding_box(lx / 2 - EPS, -EPS, -EPS, lx + EPS, +EPS, +EPS, 1)
top_right = gmsh.model.occ.get_entities_in_bounding_box(lx / 2 - EPS, ly - EPS, -EPS, lx + EPS, ly + EPS, +EPS, 1)

gmsh.model.mesh.set_periodic(1, [bot_left[1][2], bot_right[1][2]], [top_left[1][2], top_right[1][2]], y_translation[:])

left_bot = gmsh.model.occ.get_entities_in_bounding_box(-EPS, -EPS, -EPS, + EPS, ly/2+EPS, +EPS, 1)
left_top = gmsh.model.occ.get_entities_in_bounding_box(-EPS, ly/2 - EPS, -EPS, + EPS, ly + EPS, +EPS, 1)
right_bot = gmsh.model.occ.get_entities_in_bounding_box(lx - EPS, -EPS, -EPS, lx + EPS, ly/2+EPS, +EPS, 1)
right_top = gmsh.model.occ.get_entities_in_bounding_box(lx - EPS, ly/2 - EPS, -EPS, lx + EPS, ly + EPS, +EPS, 1)

gmsh.model.mesh.set_periodic(1, [left_bot[1][2], left_top[1][2]], [right_bot[1][2], right_top[1][2]], x_translation[:])

# Физические группы (границы, углы и объем)
left_central_nodes = gmsh.model.occ.get_entities_in_bounding_box(-EPS, +EPS, -EPS, +EPS, ly - EPS, +EPS, 0)
right_central_nodes = gmsh.model.occ.get_entities_in_bounding_box(lx-EPS, +EPS, -EPS, lx+EPS, ly - EPS, +EPS, 0)
bottom_central_nodes = gmsh.model.occ.get_entities_in_bounding_box(+EPS, -EPS, -EPS, lx-EPS, +EPS, +EPS, 0)
top_central_nodes = gmsh.model.occ.get_entities_in_bounding_box(+EPS, ly-EPS, -EPS, lx-EPS, ly+EPS, +EPS, 0)

left_bottom_node = gmsh.model.occ.get_entities_in_bounding_box(-EPS, -EPS, -EPS, +EPS, +EPS, +EPS, 0)
right_bottom_node = gmsh.model.occ.get_entities_in_bounding_box(lx-EPS, -EPS, -EPS, lx+EPS, +EPS, +EPS, 0)
left_top_node = gmsh.model.occ.get_entities_in_bounding_box(-EPS, ly-EPS, -EPS, +EPS, ly+EPS, +EPS, 0)
right_top_node = gmsh.model.occ.get_entities_in_bounding_box(lx-EPS, ly-EPS, -EPS, lx+EPS, ly+EPS, +EPS, 0)

gmsh.model.add_physical_group(1, [left_bot[1][2], left_top[1][2]], -1, "left")
gmsh.model.add_physical_group(0, [dimtag[2] for dimtag in left_central_nodes], -1, "left")
gmsh.model.add_physical_group(1, [right_bot[1][2], right_top[1][2]], -1, "right")
gmsh.model.add_physical_group(0, [dimtag[2] for dimtag in right_central_nodes], -1, "right")
gmsh.model.add_physical_group(1, [bot_left[1][2], bot_right[1][2]], -1, "bottom")
gmsh.model.add_physical_group(0, [dimtag[2] for dimtag in bottom_central_nodes], -1, "bottom")
gmsh.model.add_physical_group(1, [top_left[1][2], top_right[1][2]], -1, "top")
gmsh.model.add_physical_group(0, [dimtag[2] for dimtag in top_central_nodes], -1, "top")

gmsh.model.add_physical_group(0, [left_top_node[1][2]], -1, "left_top")
gmsh.model.add_physical_group(0, [right_top_node[1][2]], -1, "right_top")
gmsh.model.add_physical_group(0, [left_bottom_node[1][2]], -1, "left_bottom")
gmsh.model.add_physical_group(0, [right_bottom_node[1][2]], -1, "right_bottom")

main_body = gmsh.model.occ.get_entities(2)
gmsh.model.add_physical_group(2, [main_body[1][2]], -1, "rve")

gmsh.model.mesh.setSize(gmsh.model.getEntities(0), t * 0.1)
gmsh.option.setNumber("Mesh.Algorithm", 6)
gmsh.model.mesh.generate(2)

mesh_file = (@__DIR__) * "/honeycomb_rve.msh"
gmsh.write(mesh_file)
gmsh.finalize()

Здесь одним из ключевых шагов является использование функции set_periodic в Gmsh, которая гарантирует, что узлы сетки на левой границе будут иметь точные зеркальные копии на правой (аналогично для верхней и нижней границ). Для определения периодичности необходимо содать матрицу трансформации $\mathbf{T} \in \mathbb{R}^{4 \times 4}$ в гомогенных координатах (x_translation и y_translation для пар лево-право и верх-низ соответственно), которая сочетает в себе поворот и перемещение:

\[ \mathbf{T}=\begin{bmatrix} \mathbf{R}&\mathbf{v}\\ \mathbf{0} & 1 \end{bmatrix}, \]

где $\mathbf{R}$ - матрица поворота (см. Глава 4) и $\mathbf{v}$ – вектор перемещений.

5.3.3 Подготовка пространств и базовой матрицы

Поскольку мы будем использовать метод множителей Лагранжа для наложения граничных условий, то мы не накладываем условия Дирихле на пространствах Gridap. Вместо этого, мы собираем полную “свободную” матрицу жесткости, а затем явно извлекаем узлы для связывания. Чтобы связать $i$-й левый узел с $i$-м правым, списки узлов нужно отсортировать по координатам: для левой и правой грани необходимо отсортировать по $y$ координате, а для верней и нижней - по $x$.

Список 5.2: Сборка свободной матрицы жесткости и сортировка узлов

model = GmshDiscreteModel(mesh_file)

left_nodes = get_nodes_by_tag(model, "left")
right_nodes = get_nodes_by_tag(model, "right")
top_nodes = get_nodes_by_tag(model, "top")
bottom_nodes = get_nodes_by_tag(model, "bottom")

lb_node = get_nodes_by_tag(model, "left_bottom")
rb_node = get_nodes_by_tag(model, "right_bottom")
lt_node = get_nodes_by_tag(model, "left_top")
rt_node = get_nodes_by_tag(model, "right_top")

# Сортировка узлов для корректного наложения периодичности
function reorder!(nodes, dir, coordinates)
    projections = [crd[dir] for crd in coordinates[nodes]]
    nodes .= nodes[sortperm(projections)]
end

reorder!(left_nodes, 2, model.grid.node_coordinates)
reorder!(right_nodes, 2, model.grid.node_coordinates)
reorder!(top_nodes, 1, model.grid.node_coordinates)
reorder!(bottom_nodes, 1, model.grid.node_coordinates)

E = 70e9
ν = 0.33
λ = (E * ν) / ((1 + ν) * (1 - 2 * ν))
μ = E / (2 * (1 + ν))
σ(ε) = λ * tr(ε) * one(ε) + 2 * μ * ε

reffe = ReferenceFE(lagrangian, VectorValue{2,Float64}, 1)
V = TestFESpace(model, reffe, conformity=:H1)
U = TrialFESpace(V) # Без ГУ

EL_ORDER = 1
degree = EL_ORDER * 2
Ω = Triangulation(model)
dΩ = Measure(Ω, degree)

blf_stiff(u, v) = ∫(ε(v) ⊙ (σ ∘ ε(u))) * dΩ
mK = assemble_matrix(blf_stiff, V, U)

5.3.4 Формирование матрицы ограничений (PBC)

Мы создаем разреженную матрицу $\mathbf{A}$, в которой каждая строка отвечает за одно уравнение связи. Например, уравнение $\mathbf{u}_{\text{right}} - \mathbf{u}_{\text{left}} = \bar{\varepsilon}_{xx} \cdot l_x$ кодируется как $1 \cdot u_{rx} - 1 \cdot u_{lx} = q$.

Список 5.3: Построение матрицы связей А

num_l_to_r_dofs = 2 * length(left_nodes)
num_t_to_b_dofs = 2 * length(top_nodes)
num_corner_dofs = 2 * 4

num_constrains = num_l_to_r_dofs + num_t_to_b_dofs + num_corner_dofs
num_glob_dofs = size(mK, 1)

# Вспомогательная функция для получения DOF-индексов узлов
function get_dofs(nodes)
    dx = [V.metadata.node_and_comp_to_dof[n][1] for n in nodes]
    dy = [V.metadata.node_and_comp_to_dof[n][2] for n in nodes]
    return dx, dy
end

l_dx, l_dy = get_dofs(left_nodes); r_dx, r_dy = get_dofs(right_nodes)
b_dx, b_dy = get_dofs(bottom_nodes); t_dx, t_dy = get_dofs(top_nodes)
lb_dx, lb_dy = get_dofs(lb_node); rb_dx, rb_dy = get_dofs(rb_node)
lt_dx, lt_dy = get_dofs(lt_node); rt_dx, rt_dy = get_dofs(rt_node)

A = spzeros(num_constrains, num_glob_dofs)

# Левая-Правая и Верх-Низ периодичность ( U_right - U_left = dU )
A[1:(num_l_to_r_dofs+num_t_to_b_dofs), [l_dx; l_dy; t_dx; t_dy]] = I(num_l_to_r_dofs+num_t_to_b_dofs)
A[1:(num_l_to_r_dofs+num_t_to_b_dofs), [r_dx; r_dy; b_dx; b_dy]] = -I(num_l_to_r_dofs+num_t_to_b_dofs)

# Углы ячейки задают макроскопическую деформацию, их перемещения фиксируются явно
A[(num_l_to_r_dofs+num_t_to_b_dofs+1):end, [lb_dx; lb_dy; rb_dx; rb_dy; lt_dx; lt_dy; rt_dx; rt_dy]] = I(num_corner_dofs)

# Разметка строк матрицы A
l_to_r_x_eqs = 1:length(l_dx)
l_to_r_y_eqs = (1:length(l_dy)) .+ length(l_dx)
t_to_b_x_eqs = (1:length(t_dx)) .+ length(l_dy) .+ length(l_dx)
t_to_b_y_eqs = (1:length(t_dy)) .+ length(t_dx) .+ length(l_dy) .+ length(l_dx)
lb_x_eq, lb_y_eq, rb_x_eq, rb_y_eq, lt_x_eq, lt_y_eq, rt_x_eq, rt_y_eq = (t_to_b_y_eqs[end]+1):num_constrains

5.3.5 Решение

Решим блочную седловую задачу для 3-х вариантов нагрузки: макро-растяжение по X, по Y и макро-сдвиг.

Список 5.4: Решение седловой задачи для 3-х состояний и интеграция напряжений

# Формирование блочной матрицы M
M = [mK A'; A spzeros(num_constrains, num_constrains)]

## СЛУЧАЙ 1: Растяжение вдоль X (Деформация = 10%)
rhs = zeros(num_constrains)
rhs[l_to_r_x_eqs] .= -0.1lx; rhs[l_to_r_y_eqs] .= 0.0
rhs[t_to_b_x_eqs] .= 0.0;    rhs[t_to_b_y_eqs] .= 0.0
rhs[lb_x_eq] = 0.0; rhs[lb_y_eq] = 0.0
rhs[rb_x_eq] = 0.1lx; rhs[rb_y_eq] = 0.0
rhs[lt_x_eq] = 0.0; rhs[lt_y_eq] = 0.0
rhs[rt_x_eq] = 0.1lx; rhs[rt_y_eq] = 0.0

F = [zeros(num_glob_dofs); rhs]
sol = M \ F
uh1 = FEFunction(V, sol[1:num_glob_dofs])

## СЛУЧАЙ 2: Растяжение вдоль Y
rhs = zeros(num_constrains)
rhs[l_to_r_x_eqs] .= 0.0; rhs[l_to_r_y_eqs] .= 0.0
rhs[t_to_b_x_eqs] .= 0.0; rhs[t_to_b_y_eqs] .= 0.1*ly
rhs[lb_x_eq] = 0.0; rhs[lb_y_eq] = 0.0
rhs[rb_x_eq] = 0.0; rhs[rb_y_eq] = 0.0
rhs[lt_x_eq] = 0.0; rhs[lt_y_eq] = 0.1ly
rhs[rt_x_eq] = 0.0; rhs[rt_y_eq] = 0.1ly

F = [zeros(num_glob_dofs); rhs]
sol = M \ F
uh2 = FEFunction(V, sol[1:num_glob_dofs])

## СЛУЧАЙ 3: Сдвиг в плоскости XY
rhs = zeros(num_constrains)
rhs[l_to_r_x_eqs] .= 0.0; rhs[l_to_r_y_eqs] .= -0.1lx
rhs[t_to_b_x_eqs] .= 0.1ly; rhs[t_to_b_y_eqs] .= 0.0
rhs[lb_x_eq] = 0.0; rhs[lb_y_eq] = 0.0
rhs[rb_x_eq] = 0.0; rhs[rb_y_eq] = 0.1lx
rhs[lt_x_eq] = 0.1ly; rhs[lt_y_eq] = 0.0
rhs[rt_x_eq] = 0.1ly; rhs[rt_y_eq] = 0.1lx

F = [zeros(num_glob_dofs); rhs]
sol = M \ F
uh3 = FEFunction(V, sol[1:num_glob_dofs])

# Осреднение напряжений
σ₁ = sum(∫(σ ∘ ε(uh1)) * dΩ)
σ₂ = sum(∫(σ ∘ ε(uh2)) * dΩ)
σ₃ = sum(∫(σ ∘ ε(uh3)) * dΩ)

Gridap.TensorValues.SymTensorValue{2, Float64, 3}(-4957.554681112131, 1.2658558852256525e9, 66234.5173207121)

5.3.6 Визуализация и интерпретация результатов

Для визуализации мы будем использовать напряжения по Мизесу. Для этого объявим функцию mses, которая вычисляет эквивалентные напряжения по мизесу из полного тензора напряжений. Далее

Рисунок 5.2: деформированное состояние ячейки и эквивалентные напряжения

Ключевые физические наблюдения:

Искривление границ (PBC): В отличие от жестких кинематических условий (KUBC), где границы ячейки оставались бы идеально прямыми линиями, при наложении PBC внешние границы RVE изгибаются. При этом искривление левой границы в точности повторяет профиль правой. Это позволяет ячейкам собираться в бесконечную “мозаику” без образования зазоров или нахлестов.

Снятие краевых эффектов: Свободное искривление границ позволяет материалу занять наиболее энергетически выгодное положение. Именно поэтому PBC обеспечивает самую точную оценку эффективной жесткости структуры без ее искусственного завышения (как в KUBC).

5.4 Контрольные вопросы и задания

5.4.1 Теоретические вопросы

Почему кинематические однородные граничные условия (KUBC) завышают эффективную жесткость материала, а статические (SUBC) — занижают?
В чем математический смысл использования метода множителей Лагранжа и блочной матрицы при наложении PBC? Как формируется матрица $\mathbf{A}$?
Какие требования предъявляются к геометрии и сетке представительного элемента объема (RVE) для корректного применения PBC?

5.4.2 Практические задания

Измените параметр пористости базовой ячейки (толщину стенок t и радиус отверстия R_hole) в скрипте генерации геометрии. Постройте график зависимости эффективной компоненты упругости $C_{11}$ (основываясь на векторе $\sigma_1$) от относительной плотности ячейки.
Исследование границ Хилла. Замените метод наложения граничных условий в коде с периодического (PBC) на кинематическое (KUBC) — для этого можно использовать стандартный механизм dirichlet_tags в Gridap, фиксируя узлы. Сравните полученные значения $\mathbf{C}_{\mathrm{eff}}$. Убедитесь на практике, что диагональные элементы матрицы при KUBC строго больше, чем при PBC.
Постройте зависимость упругих компонент от размера репрезентативного объемного элемента для $N=[1{,}2{,}4{,}8]$.

```{julia} #| echo: false #| output: false tcf(filename::String) = (@__DIR__) * "\\" * filename; import Pkg; Pkg.activate((@__DIR__) * "\\..\\.."); using CCMechBook using Gmsh using Gridap using GridapGmsh using LinearAlgebra using SparseArrays using GLMakie ``` # Вычисление эффективных свойств (гомогенизация) периодических структур ## Аннотация и цели обучения - **Предварительные знания**: Основы линейной теории упругости, нотация Фойгта (Урок 2), навыки работы с глобальными матрицами жесткости, тензорная алгебра. - **Освоенные численные методы**: Метод асимптотической гомогенизации, концепция представительного элемента объема (RVE), генерация строго периодических сеток в Gmsh, реализация периодических граничных условий (PBC) методом множителей Лагранжа через матрицу алгебраических ограничений. - **Физическая задача**: Определение макроскопических эффективных упругих характеристик (тензора упругости) гетерогенной периодической сотовой структуры, применяемой в аддитивных технологиях и при создании метаматериалов. ## Осовные теоретические сведения ### Асимптотическая гомогенизация и масштабный переход Многие современные материалы, такие как композиты, пенопласты, 3D-печатные сотовые структуры, имеют ярко выраженную внутреннюю микроструктуру. Моделировать каждую пору или каждое волокно в масштабах крупной детали (например, крыла самолета) вычислительно невозможно. Метод асимптотической гомогенизации решает эту проблему путем разделения задачи на два масштаба: 1. Микромасштаб ($l$): Уровень отдельной ячейки материала, где явно моделируются все пустоты, включения и геометрические особенности. На этом уровне материалы компонент считаются сплошными с известными свойствами. 2. Макромасштаб ($L$): Уровень всей конструкции ($l \ll L$). На этом уровне материал считается однородным (гомогенным), но обладающим эффективными (усредненными) свойствами, вычисленными из микромасштаба. Связующим звеном выступает **представительный объемный элемент** (Representative Volume Element, RVE). Для структур с периодической геометрией (например, сотовых) RVE представляет собой одну минимальную повторяющуюся ячейку. Фундаментальным условием корректного перехода между масштабами является теорема Хилла - Манделя, которая гласит, что макроскопическая работа деформации должна быть равна среднему по объему значению микроскопической работы: $$ \langle \mathbf{\sigma} : \mathbf{\varepsilon} \rangle = \langle \mathbf{\sigma} \rangle : \langle \mathbf{\varepsilon} \rangle $$ где скобки $\langle \cdot \rangle$ обозначают осреднение величины по объему RVE: $\langle f \rangle = \frac{1}{V}\int_V f \, dV$. ### Осреднение и краевые задачи На микроуровне внутри RVE (область $\Omega$) справедливы классические уравнения линейной теории упругости: $$ \nabla \cdot \mathbf{\sigma} = 0, \quad \mathbf{\sigma} = {}^4\mathbf{C} : \mathbf{\varepsilon}(\mathbf{u}) $$ После процедуры осреднения, связь между макроскопическими (эффективными) напряжениями и деформациями принимает вид обобщенного закона Гука: $$ \langle \mathbf{\sigma} \rangle = {}^4\mathbf{C}^{\mathrm{eff}} \langle \mathbf{\varepsilon} \rangle $$ где ${}^4\mathbf{C}^{\mathrm{eff}}$ - тензор эффективных упругих компонент четвертого порядка. Однако часто его представляют в виде матрицы используя нотацию Фойгта (см. Урок 2) и в двумерном случае матрица эффективных компонент $\mathbf{C}_v^{\mathrm{eff}}$ имеет размерность $3 \times 3$. Чтобы найти все её компоненты, нам необходимо решить три независимые краевые задачи, прикладывая к ячейке три макроскопических состояния деформации $\bar{\varepsilon}$: 1. Одноосное растяжение вдоль оси X: $\bar{\varepsilon}_{11} = 1$ 2. Одноосное растяжение вдоль оси Y: $\bar{\varepsilon}_{22} = 1$ 3. Чистый сдвиг в плоскости XY: $\bar{\varepsilon}_{12} = 1$ Прикладывая единичную макро-деформацию, осредненный вектор напряжений $\langle \mathbf{\sigma}_v \rangle$ будет в точности равен соответствующему столбцу матрицы $\mathbf{C}_v^{\mathrm{eff}}$. ### Типы граничных условий Существует три способа приложить макро-деформацию $\bar{\varepsilon}$ к границе ячейки: 1. Кинематические (KUBC): Жесткое задание перемещений на границе. Завышает жесткость. 2. Статические (SUBC): Задание усилий (тракций) на границе. Занижает жесткость. 3. Периодические (PBC): Противоположные границы ячейки могут искривляться свободно, но их перемещения строго связаны друг с другом: $$ \mathbf{u}(\mathbf{x}^+) - \mathbf{u}(\mathbf{x}^-) = \bar{\varepsilon} \cdot (\mathbf{x}^+ - \mathbf{x}^-) $$ где $\mathbf{x}^+$ и $\mathbf{x}^-$ — координаты соответствующих точек на противоположных гранях. Для строго периодических сред метод PBC дает точное значение эффективной жесткости. Математически доказано, что $\mathbf{C}^{SUBC} \le \mathbf{C}^{PBC} \le \mathbf{C}^{KUBC}$. Для реализации PBC нам необходимо построить идеальную периодическую геометрию и связать степени свободы противоположных узлов. Это приведет к системе уравнений с ограничениями (алгебраическими связями). В методе конечных элементов такие задачи решаются с помощью множителей Лагранжа, где исходная матрица жесткости $\mathbf{K}$ дополняется матрицей ограничений $\mathbf{A}$: $$ \begin{bmatrix} \mathbf{K} & \mathbf{A}^T \\ \mathbf{A} & \mathbf{0} \end{bmatrix} \begin{Bmatrix} \mathbf{u} \\ \mathbf{\lambda} \end{Bmatrix}= \begin{Bmatrix} \mathbf{0} \\ \mathbf{q} \end{Bmatrix} $$ Вектор $\mathbf{q}$ задает относительные смещения границ, определяемые макроскопической деформацией $\bar{\varepsilon}$. ## Практический задача ### Исходные данные Смоделируем 2D-ячейку полимерного материала полученную при помощи 3D-печати (сотовая структура). Материал: $E =$ 70 ГПа, $\nu =$ 0.33. ### Построение геометрии и периодической сетки в Gmsh Репрезентативная ячейка сотовой структуры (@fig-honeycomb) представляет собой прямоугольник с вырезанными половинками правильных шестиугольников по краям, в зависимости от того, как определить границы репрезентативной ячейки (см. @fig-rve). ::: {#fig-geometry layout-ncol=2} ![](../../assets/pics/honeycomb.png){#fig-honeycomb} ![](../../assets/pics/rve.png){#fig-rve} Репрезентативная ячейка сотовой структуры и ее параметры ::: Для построения геометрии ячейки, создадим вспомогательную функцию `create_hexagon`, которая будет создавать шестиугольник. Далее создаем прямоугольник $l_x \times l_y$ и 4 шестиугольника, которые в дальнейшем вырезаются из прямоугольника помощи булевой операции `cut` (см. @lst-gmsh-rve). ```{julia} #| echo: true #| output: false #| lst-cap: "Генерация геометрии Honeycomb и наложение периодичности сетки" #| lst-label: lst-gmsh-rve gmsh.initialize() gmsh.option.setNumber("General.Terminal", 1) gmsh.model.add("Honeycomb_RVE") R_pitch = 1.0 # Шаг сот t = 0.3 # Толщина стенки напечатанного материала R_hole = R_pitch - t / sqrt(3) lx = sqrt(3) * R_pitch # Ширина RVE ly = 3.0 * R_pitch # Высота RVE # Функция создания шестиугольника function create_hexagon(x, y, R) pts = [] # Поворачиваем на 30 градусов, чтобы грани были горизонтальными for i in 0:5 angle = i * pi / 3 + pi / 6 px = x + R * cos(angle) py = y + R * sin(angle) push!(pts, gmsh.model.occ.addPoint(px, py, 0)) end lines = [] for i in 1:6 p1 = pts[i] p2 = pts[i == 6 ? 1 : i + 1] push!(lines, gmsh.model.occ.addLine(p1, p2)) end cl = gmsh.model.occ.addCurveLoop(lines) return gmsh.model.occ.addPlaneSurface([cl]) end rect = gmsh.model.occ.addRectangle(0, 0, 0, lx, ly) hex_tags = [ create_hexagon(lx / 2, 0, R_hole), create_hexagon(lx / 2, ly, R_hole), create_hexagon(0, ly / 2, R_hole), create_hexagon(lx, ly / 2, R_hole), ] gmsh.model.occ.cut([(2, rect)], [(2, tag) for tag in hex_tags]) gmsh.model.occ.synchronize() EPS = 1e-6 # Матрицы трансляции для периодической сетки x_translation = [1 0 0 lx; 0 1 0 0; 0 0 1 0; 0 0 0 1] y_translation = [1 0 0 0; 0 1 0 ly; 0 0 1 0; 0 0 0 1] # Поиск парных границ bot_left = gmsh.model.occ.get_entities_in_bounding_box(-EPS, -EPS, -EPS, lx / 2 + EPS, +EPS, +EPS, 1) top_left = gmsh.model.occ.get_entities_in_bounding_box(-EPS, ly - EPS, -EPS, lx / 2 + EPS, ly + EPS, +EPS, 1) bot_right = gmsh.model.occ.get_entities_in_bounding_box(lx / 2 - EPS, -EPS, -EPS, lx + EPS, +EPS, +EPS, 1) top_right = gmsh.model.occ.get_entities_in_bounding_box(lx / 2 - EPS, ly - EPS, -EPS, lx + EPS, ly + EPS, +EPS, 1) gmsh.model.mesh.set_periodic(1, [bot_left[1][2], bot_right[1][2]], [top_left[1][2], top_right[1][2]], y_translation[:]) # <1> left_bot = gmsh.model.occ.get_entities_in_bounding_box(-EPS, -EPS, -EPS, + EPS, ly/2+EPS, +EPS, 1) left_top = gmsh.model.occ.get_entities_in_bounding_box(-EPS, ly/2 - EPS, -EPS, + EPS, ly + EPS, +EPS, 1) right_bot = gmsh.model.occ.get_entities_in_bounding_box(lx - EPS, -EPS, -EPS, lx + EPS, ly/2+EPS, +EPS, 1) right_top = gmsh.model.occ.get_entities_in_bounding_box(lx - EPS, ly/2 - EPS, -EPS, lx + EPS, ly + EPS, +EPS, 1) gmsh.model.mesh.set_periodic(1, [left_bot[1][2], left_top[1][2]], [right_bot[1][2], right_top[1][2]], x_translation[:]) # <2> # Физические группы (границы, углы и объем) left_central_nodes = gmsh.model.occ.get_entities_in_bounding_box(-EPS, +EPS, -EPS, +EPS, ly - EPS, +EPS, 0) right_central_nodes = gmsh.model.occ.get_entities_in_bounding_box(lx-EPS, +EPS, -EPS, lx+EPS, ly - EPS, +EPS, 0) bottom_central_nodes = gmsh.model.occ.get_entities_in_bounding_box(+EPS, -EPS, -EPS, lx-EPS, +EPS, +EPS, 0) top_central_nodes = gmsh.model.occ.get_entities_in_bounding_box(+EPS, ly-EPS, -EPS, lx-EPS, ly+EPS, +EPS, 0) left_bottom_node = gmsh.model.occ.get_entities_in_bounding_box(-EPS, -EPS, -EPS, +EPS, +EPS, +EPS, 0) right_bottom_node = gmsh.model.occ.get_entities_in_bounding_box(lx-EPS, -EPS, -EPS, lx+EPS, +EPS, +EPS, 0) left_top_node = gmsh.model.occ.get_entities_in_bounding_box(-EPS, ly-EPS, -EPS, +EPS, ly+EPS, +EPS, 0) right_top_node = gmsh.model.occ.get_entities_in_bounding_box(lx-EPS, ly-EPS, -EPS, lx+EPS, ly+EPS, +EPS, 0) gmsh.model.add_physical_group(1, [left_bot[1][2], left_top[1][2]], -1, "left") gmsh.model.add_physical_group(0, [dimtag[2] for dimtag in left_central_nodes], -1, "left") gmsh.model.add_physical_group(1, [right_bot[1][2], right_top[1][2]], -1, "right") gmsh.model.add_physical_group(0, [dimtag[2] for dimtag in right_central_nodes], -1, "right") gmsh.model.add_physical_group(1, [bot_left[1][2], bot_right[1][2]], -1, "bottom") gmsh.model.add_physical_group(0, [dimtag[2] for dimtag in bottom_central_nodes], -1, "bottom") gmsh.model.add_physical_group(1, [top_left[1][2], top_right[1][2]], -1, "top") gmsh.model.add_physical_group(0, [dimtag[2] for dimtag in top_central_nodes], -1, "top") gmsh.model.add_physical_group(0, [left_top_node[1][2]], -1, "left_top") gmsh.model.add_physical_group(0, [right_top_node[1][2]], -1, "right_top") gmsh.model.add_physical_group(0, [left_bottom_node[1][2]], -1, "left_bottom") gmsh.model.add_physical_group(0, [right_bottom_node[1][2]], -1, "right_bottom") main_body = gmsh.model.occ.get_entities(2) gmsh.model.add_physical_group(2, [main_body[1][2]], -1, "rve") gmsh.model.mesh.setSize(gmsh.model.getEntities(0), t * 0.1) gmsh.option.setNumber("Mesh.Algorithm", 6) gmsh.model.mesh.generate(2) mesh_file = (@__DIR__) * "/honeycomb_rve.msh" gmsh.write(mesh_file) gmsh.finalize() ``` Здесь одним из ключевых шагов является использование функции `set_periodic` в **Gmsh**, которая гарантирует, что узлы сетки на левой границе будут иметь точные зеркальные копии на правой (аналогично для верхней и нижней границ). Для определения периодичности необходимо содать матрицу трансформации $\mathbf{T} \in \mathbb{R}^{4 \times 4}$ в гомогенных координатах (`x_translation` и `y_translation` для пар *лево-право* и *верх-низ* соответственно), которая сочетает в себе поворот и перемещение: $$ \mathbf{T}=\begin{bmatrix} \mathbf{R}&\mathbf{v}\\ \mathbf{0} & 1 \end{bmatrix}, $$ где $\mathbf{R}$ - матрица поворота (см. @sec-modal-analysis) и $\mathbf{v}$ -- вектор перемещений. ### Подготовка пространств и базовой матрицы Поскольку мы будем использовать метод множителей Лагранжа для наложения граничных условий, то мы не накладываем условия Дирихле на пространствах Gridap. Вместо этого, мы собираем полную "свободную" матрицу жесткости, а затем явно извлекаем узлы для связывания. Чтобы связать $i$-й левый узел с $i$-м правым, списки узлов нужно отсортировать по координатам: для левой и правой грани необходимо отсортировать по $y$ координате, а для верней и нижней - по $x$. ```{julia} #| echo: true #| output: false #| lst-cap: "Сборка свободной матрицы жесткости и сортировка узлов" #| lst-label: lst-gridap-setup model = GmshDiscreteModel(mesh_file) left_nodes = get_nodes_by_tag(model, "left") right_nodes = get_nodes_by_tag(model, "right") top_nodes = get_nodes_by_tag(model, "top") bottom_nodes = get_nodes_by_tag(model, "bottom") lb_node = get_nodes_by_tag(model, "left_bottom") rb_node = get_nodes_by_tag(model, "right_bottom") lt_node = get_nodes_by_tag(model, "left_top") rt_node = get_nodes_by_tag(model, "right_top") # Сортировка узлов для корректного наложения периодичности function reorder!(nodes, dir, coordinates) projections = [crd[dir] for crd in coordinates[nodes]] nodes .= nodes[sortperm(projections)] end reorder!(left_nodes, 2, model.grid.node_coordinates) reorder!(right_nodes, 2, model.grid.node_coordinates) reorder!(top_nodes, 1, model.grid.node_coordinates) reorder!(bottom_nodes, 1, model.grid.node_coordinates) E = 70e9 ν = 0.33 λ = (E * ν) / ((1 + ν) * (1 - 2 * ν)) μ = E / (2 * (1 + ν)) σ(ε) = λ * tr(ε) * one(ε) + 2 * μ * ε reffe = ReferenceFE(lagrangian, VectorValue{2,Float64}, 1) V = TestFESpace(model, reffe, conformity=:H1) U = TrialFESpace(V) # Без ГУ EL_ORDER = 1 degree = EL_ORDER * 2 Ω = Triangulation(model) dΩ = Measure(Ω, degree) blf_stiff(u, v) = ∫(ε(v) ⊙ (σ ∘ ε(u))) * dΩ mK = assemble_matrix(blf_stiff, V, U) ``` ### Формирование матрицы ограничений (PBC) Мы создаем разреженную матрицу $\mathbf{A}$, в которой каждая строка отвечает за одно уравнение связи. Например, уравнение $\mathbf{u}_{\text{right}} - \mathbf{u}_{\text{left}} = \bar{\varepsilon}_{xx} \cdot l_x$ кодируется как $1 \cdot u_{rx} - 1 \cdot u_{lx} = q$. ```{julia} #| echo: true #| output: false #| lst-cap: "Построение матрицы связей А" #| lst-label: lst-constraint-matrix num_l_to_r_dofs = 2 * length(left_nodes) num_t_to_b_dofs = 2 * length(top_nodes) num_corner_dofs = 2 * 4 num_constrains = num_l_to_r_dofs + num_t_to_b_dofs + num_corner_dofs num_glob_dofs = size(mK, 1) # Вспомогательная функция для получения DOF-индексов узлов function get_dofs(nodes) dx = [V.metadata.node_and_comp_to_dof[n][1] for n in nodes] dy = [V.metadata.node_and_comp_to_dof[n][2] for n in nodes] return dx, dy end l_dx, l_dy = get_dofs(left_nodes); r_dx, r_dy = get_dofs(right_nodes) b_dx, b_dy = get_dofs(bottom_nodes); t_dx, t_dy = get_dofs(top_nodes) lb_dx, lb_dy = get_dofs(lb_node); rb_dx, rb_dy = get_dofs(rb_node) lt_dx, lt_dy = get_dofs(lt_node); rt_dx, rt_dy = get_dofs(rt_node) A = spzeros(num_constrains, num_glob_dofs) # Левая-Правая и Верх-Низ периодичность ( U_right - U_left = dU ) A[1:(num_l_to_r_dofs+num_t_to_b_dofs), [l_dx; l_dy; t_dx; t_dy]] = I(num_l_to_r_dofs+num_t_to_b_dofs) A[1:(num_l_to_r_dofs+num_t_to_b_dofs), [r_dx; r_dy; b_dx; b_dy]] = -I(num_l_to_r_dofs+num_t_to_b_dofs) # Углы ячейки задают макроскопическую деформацию, их перемещения фиксируются явно A[(num_l_to_r_dofs+num_t_to_b_dofs+1):end, [lb_dx; lb_dy; rb_dx; rb_dy; lt_dx; lt_dy; rt_dx; rt_dy]] = I(num_corner_dofs) # Разметка строк матрицы A l_to_r_x_eqs = 1:length(l_dx) l_to_r_y_eqs = (1:length(l_dy)) .+ length(l_dx) t_to_b_x_eqs = (1:length(t_dx)) .+ length(l_dy) .+ length(l_dx) t_to_b_y_eqs = (1:length(t_dy)) .+ length(t_dx) .+ length(l_dy) .+ length(l_dx) lb_x_eq, lb_y_eq, rb_x_eq, rb_y_eq, lt_x_eq, lt_y_eq, rt_x_eq, rt_y_eq = (t_to_b_y_eqs[end]+1):num_constrains ``` ### Решение Решим блочную седловую задачу для 3-х вариантов нагрузки: макро-растяжение по X, по Y и макро-сдвиг. ```{julia} #| echo: true #| output: true #| lst-cap: "Решение седловой задачи для 3-х состояний и интеграция напряжений" #| lst-label: lst-solve-homogenization # Формирование блочной матрицы M M = [mK A'; A spzeros(num_constrains, num_constrains)] ## СЛУЧАЙ 1: Растяжение вдоль X (Деформация = 10%) rhs = zeros(num_constrains) rhs[l_to_r_x_eqs] .= -0.1lx; rhs[l_to_r_y_eqs] .= 0.0 rhs[t_to_b_x_eqs] .= 0.0; rhs[t_to_b_y_eqs] .= 0.0 rhs[lb_x_eq] = 0.0; rhs[lb_y_eq] = 0.0 rhs[rb_x_eq] = 0.1lx; rhs[rb_y_eq] = 0.0 rhs[lt_x_eq] = 0.0; rhs[lt_y_eq] = 0.0 rhs[rt_x_eq] = 0.1lx; rhs[rt_y_eq] = 0.0 F = [zeros(num_glob_dofs); rhs] sol = M \ F uh1 = FEFunction(V, sol[1:num_glob_dofs]) ## СЛУЧАЙ 2: Растяжение вдоль Y rhs = zeros(num_constrains) rhs[l_to_r_x_eqs] .= 0.0; rhs[l_to_r_y_eqs] .= 0.0 rhs[t_to_b_x_eqs] .= 0.0; rhs[t_to_b_y_eqs] .= 0.1*ly rhs[lb_x_eq] = 0.0; rhs[lb_y_eq] = 0.0 rhs[rb_x_eq] = 0.0; rhs[rb_y_eq] = 0.0 rhs[lt_x_eq] = 0.0; rhs[lt_y_eq] = 0.1ly rhs[rt_x_eq] = 0.0; rhs[rt_y_eq] = 0.1ly F = [zeros(num_glob_dofs); rhs] sol = M \ F uh2 = FEFunction(V, sol[1:num_glob_dofs]) ## СЛУЧАЙ 3: Сдвиг в плоскости XY rhs = zeros(num_constrains) rhs[l_to_r_x_eqs] .= 0.0; rhs[l_to_r_y_eqs] .= -0.1lx rhs[t_to_b_x_eqs] .= 0.1ly; rhs[t_to_b_y_eqs] .= 0.0 rhs[lb_x_eq] = 0.0; rhs[lb_y_eq] = 0.0 rhs[rb_x_eq] = 0.0; rhs[rb_y_eq] = 0.1lx rhs[lt_x_eq] = 0.1ly; rhs[lt_y_eq] = 0.0 rhs[rt_x_eq] = 0.1ly; rhs[rt_y_eq] = 0.1lx F = [zeros(num_glob_dofs); rhs] sol = M \ F uh3 = FEFunction(V, sol[1:num_glob_dofs]) # Осреднение напряжений σ₁ = sum(∫(σ ∘ ε(uh1)) * dΩ) σ₂ = sum(∫(σ ∘ ε(uh2)) * dΩ) σ₃ = sum(∫(σ ∘ ε(uh3)) * dΩ) ``` ### Визуализация и интерпретация результатов Для визуализации мы будем использовать напряжения по Мизесу. Для этого объявим функцию `mses`, которая вычисляет эквивалентные напряжения по мизесу из полного тензора напряжений. Далее ```{julia} #| echo: false #| output: true #| label: fig-result #| fig-cap: "деформированное состояние ячейки и эквивалентные напряжения" function mses(σ) d = size(σ)[1] I = one(σ) s = σ - (1/3) * tr(σ) * I return sqrt(3/2 * inner(s, s)) end fig = Figure(size=(700, 360)) ax1 = Axis(fig[1, 1], aspect=DataAspect(), title="а)", titlealign=:left) plt1 = warpedmeshtri!(ax1, Ω, warp=uh1, colorby=(x)->(mses ∘ σ ∘ ε(uh1))(x), colormap=:viridis, nomesh=true) ax2 = Axis(fig[1, 2], aspect=DataAspect(), title="b)", titlealign=:left) plt2 = warpedmeshtri!(ax2, Ω, warp=uh2, colorby=(x)->(mses ∘ σ ∘ ε(uh2))(x), colormap=:viridis, nomesh=true) ax3 = Axis(fig[1, 3], aspect=DataAspect(), title="c)", titlealign=:left) plt3 = warpedmeshtri!(ax3, Ω, warp=uh3, colorby=(x)->(mses ∘ σ ∘ ε(uh3))(x), colormap=:viridis, nomesh=true) fig ``` Ключевые физические наблюдения: **Искривление границ (PBC)**: В отличие от жестких кинематических условий (KUBC), где границы ячейки оставались бы идеально прямыми линиями, при наложении PBC внешние границы RVE изгибаются. При этом искривление левой границы в точности повторяет профиль правой. Это позволяет ячейкам собираться в бесконечную "мозаику" без образования зазоров или нахлестов. **Снятие краевых эффектов**: Свободное искривление границ позволяет материалу занять наиболее энергетически выгодное положение. Именно поэтому PBC обеспечивает самую точную оценку эффективной жесткости структуры без ее искусственного завышения (как в KUBC). ## Контрольные вопросы и задания ### Теоретические вопросы 1. Почему кинематические однородные граничные условия (KUBC) завышают эффективную жесткость материала, а статические (SUBC) — занижают? 2. В чем математический смысл использования метода множителей Лагранжа и блочной матрицы при наложении PBC? Как формируется матрица $\mathbf{A}$? 3. Какие требования предъявляются к геометрии и сетке представительного элемента объема (RVE) для корректного применения PBC? ### Практические задания 1. Измените параметр пористости базовой ячейки (толщину стенок t и радиус отверстия R_hole) в скрипте генерации геометрии. Постройте график зависимости эффективной компоненты упругости $C_{11}$ (основываясь на векторе $\sigma_1$) от относительной плотности ячейки. 2. Исследование границ Хилла. Замените метод наложения граничных условий в коде с периодического (PBC) на кинематическое (KUBC) — для этого можно использовать стандартный механизм dirichlet_tags в Gridap, фиксируя узлы. Сравните полученные значения $\mathbf{C}_{\mathrm{eff}}$. Убедитесь на практике, что диагональные элементы матрицы при KUBC строго больше, чем при PBC. 3. Постройте зависимость упругих компонент от размера репрезентативного объемного элемента для $N=[1{,}2{,}4{,}8]$.