2 Моделирование растяжения плоского образца

2.1 Аннотация

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

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

2.2.1 Упругое деформирование

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

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

\[ \begin{cases} \varepsilon_x = \frac{1}{E}\left[\sigma_x - \nu (\sigma_y + \sigma_z)\right]; \\ \varepsilon_y = \frac{1}{E}\left[\sigma_y - \nu (\sigma_x + \sigma_z)\right]; \\ \varepsilon_z = \frac{1}{E}\left[\sigma_z - \nu (\sigma_x + \sigma_y)\right]; \\ \gamma_{xy} = \frac{\tau_{xy}}{G}; \\ \gamma_{yz} = \frac{\tau_{yz}}{G}; \\ \gamma_{zx} = \frac{\tau_{zx}}{G}; \end{cases} \]

где $\varepsilon$ — линейная деформация вдоль оси; $\sigma$ — нормальное напряжение, Па; $E$ — модуль Юнга, Па; $\nu$ — коэффициент Пуассона, б.р.; $\tau$ — касательное напряжение, Па; $\gamma$ — угловая деформация, б.р.; $G$ — модуль сдвига, $G = \frac{E}{2(1+\nu)}$, Па.

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

2.2.2 Нотация Фойгта

Закон Гука выше можно представить в виде линейного отображения: \[ \sigma = {}^4\mathbf{C}_{\mathrm{eff}}{:}\varepsilon \]

где $\sigma$ и $\varepsilon$ - тензоры напряжений и деформаций, а ${}^4\mathbf{C}$ - тензор упругих постоянных. Поскольку напряжения и деформации являются тензорами второго порядка, отсюда следует что ${}^4\mathbf{C}$ является тензором четвертого порядка. Однако в силу того, что такое представление не очень удобно для записи, широкое распространение получила нотация Фойгта. Согласно этой нотации, если учесть симметрию тензоров напряжений и деформаций, то их можно представить в виде векторов, содержащих 6 компонент:

\[ \sigma_v = \begin{Bmatrix} \sigma_{11} \\ \sigma_{22} \\ \sigma_{33} \\ \sigma_{23} \\ \sigma_{13} \\ \sigma_{12} \end{Bmatrix}= \begin{Bmatrix} \sigma_x \\ \sigma_y \\ \sigma_z \\ \tau_{yz} \\ \tau_{xz} \\ \tau_{xy} \end{Bmatrix} \qquad \varepsilon_v = \begin{Bmatrix} \varepsilon_{11} \\ \varepsilon_{22} \\ \varepsilon_{33} \\ 2\varepsilon_{23} \\ 2\varepsilon_{13} \\ 2\varepsilon_{12} \end{Bmatrix}= \begin{Bmatrix} \varepsilon_x \\ \varepsilon_y \\ \varepsilon_z \\ \gamma_{yz} \\ \gamma_{xz} \\ \gamma_{xy} \end{Bmatrix} \]

Тогда тензор упругих постоянных примет вид матрицы:

\[ \mathbf{C}_v = \frac{E(1-\nu)}{(1+\nu)(1-2\nu)} \begin{bmatrix} 1 & \frac{\nu}{1-\nu} & \frac{\nu}{1-\nu} & 0 & 0 & 0 \\ \frac{\nu}{1-\nu} & 1 & \frac{\nu}{1-\nu} & 0 & 0 & 0 \\ \frac{\nu}{1-\nu} & \frac{\nu}{1-\nu} & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{1-2\nu}{2(1-\nu)} & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{1-2\nu}{2(1-\nu)} & 0 \\ 0 & 0 & 0 & 0 & 0 & \frac{1-2\nu}{2(1-\nu)} \end{bmatrix} \]

Сам закон Гука в матричном виде записывается компактно:

\[ \sigma_v = \mathbf{C}_v \varepsilon_v \]

2.2.3 Плоское напряженное и плоское деформированное состояния

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

Плоское напряженное состояние — напряжения вдоль третьей оси полагаются нулевыми, а деформации свободны (тело имеет малую толщину, например, тонкая пластина).
Плоское деформированное состояние — деформации (и перемещения) вдоль третьей оси жестко зафиксированы, а напряжения остаются свободными (тело имеет большую толщину, например, протяженная плотина).

Для плоского деформированного состояния матрица упругости имеет вид:

\[ \mathbf{C}_v = \frac{E(1-\nu)}{(1+\nu)(1-2\nu)} \begin{bmatrix} 1 & \frac{\nu}{1-\nu} & 0 \\ \frac{\nu}{1-\nu} & 1 & 0 \\ 0 & 0 & \frac{1-2\nu}{2(1-\nu)} \end{bmatrix} \]

Для плоского напряженного состояния:

\[ \mathbf{C}_{v} = \frac{E}{1-\nu^2} \begin{bmatrix} 1 & \nu & 0 \\ \nu & 1 & 0 \\ 0 & 0 & \frac{1-\nu}{2} \end{bmatrix} \]

2.2.4 Слабая формулировка задачи упругости

Рассмотрим упругое тело, занимающее область $\Omega$. В состоянии статического равновесия для него справедливо уравнение Навье:

\[ -\nabla \cdot \sigma(u) = f \quad \text{в} \quad \Omega \]

где $\sigma(u)$ — тензор напряжений, $f$ — вектор объемных сил. Для получения слабой формы умножим обе части уравнения на произвольную векторную пробную функцию $v$ и проинтегрируем по объему $\Omega$:

\[ -\int_{\Omega} (\nabla \cdot \sigma(u)) \cdot v \, dx = \int_{\Omega} f \cdot v \, dx \]

Применяя интегрирование по частям (формулу Гаусса-Остроградского) к члену с дивергенцией, получаем:

\[ -\int_{\Omega} (\nabla \cdot \sigma(u)) \cdot v \, dx = \int_{\Omega} \sigma(u) \cdot \nabla v \, dx - \int_{\partial\Omega} (\sigma(u) \cdot n) \cdot v \, dx \]

где $n$ — внешняя нормаль к границе $\partial\Omega$. Величина $\sigma(u) \cdot n$ представляет собой вектор поверхностных сил (напряжений) на границе.

Так как тензор напряжений $\sigma$ симметричен, его внутреннее произведение с антисимметричной частью градиента пробной функции $\nabla v$ равно нулю. Поэтому полный градиент $\nabla v$ можно заменить на его симметричную часть $\varepsilon(v)$ (тензор деформаций пробного поля). Это приводит к билинейной форме:

\[ a(u, v) = \int_{\Omega} \sigma(u) \cdot \varepsilon(v) \, dx \]

Линейная форма (правая часть уравнения) с учетом граничных условий Неймана $\Gamma_N$ примет вид:

\[ L(v) = \int_{\Omega} f \cdot v \, dx + \int_{\Gamma_N} \bar{t} \cdot v \, ds \]

где $\bar{t}$ — заданный вектор поверхностных нагрузок. В нашем случае мы будем рассматривать задачу без объемных сил ($f = 0$) и без распределенных поверхностных нагрузок ($\bar{t} = 0$), применяя только кинематические граничные условия (перемещения).

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

2.3.1 Постановка задачи

Смоделируем растяжение плоского стального образца. Геометрия соответствует стандартному плоскому образцу по ГОСТ 1497-84 (тип I, №6). На левой границе образца перемещения жестко зафиксированы (защемление), а правая граница смещается на 1 мм вдоль продольной оси.

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

Материал: Сталь
Модуль Юнга: E = 2.1·10¹¹ Па
Коэффициент Пуассона: ν = 0,3
Толщина образца: B = 0,04 м
Начальная длина: l₀ = 0,14 м

2.3.2 Реализация в программном стеке

2.3.2.1 Создание геометрии в Gmsh

В этом блоке мы формируем параметрическую 2D модель плоского образца со скруглениями (галтелями) и генерируем сетку.

Список 2.1

order = 1
gmsh.initialize()

b_width = 0.04
h1 = 0.08
r = (b_width - 0.03) / 2
l = 0.14 + 2 * sqrt(0.02 * 0.03)
lc = 0.005

gmsh.option.set_number("General.Terminal", 1)
gmsh.option.set_number("Mesh.ElementOrder", order)
gmsh.option.set_number("Mesh.Algorithm", 5)

gmsh.model.add("flat_sample")
factory = gmsh.model.geo

# Точки контура
factory.add_point(0, 0, 0, lc, 1)
factory.add_point(h1, 0, 0, lc, 2)
factory.add_point(h1 + r, 0, 0, lc, 3)
factory.add_point(h1 + r, r, 0, lc, 4)
factory.add_point(h1 + r + l, r, 0, lc, 5)
factory.add_point(h1 + r + l, 0, 0, lc, 6)
factory.add_point(h1 + r + l + r, 0, 0, lc, 7)
factory.add_point(h1 + r + l + r + h1, 0, 0, lc, 8)
factory.add_point(h1 + r + l + r + h1, b_width, 0, lc, 9)
factory.add_point(h1 + r + l + r, b_width, 0, lc, 10)
factory.add_point(h1 + r + l, b_width, 0, lc, 11)
factory.add_point(h1 + r + l, b_width - r, 0, lc, 12)
factory.add_point(h1 + r, b_width - r, 0, lc, 13)
factory.add_point(h1 + r, b_width, 0, lc, 14)
factory.add_point(h1, b_width, 0, lc, 15)
factory.add_point(0, b_width, 0, lc, 16)

# Линии и дуги
factory.add_line(1, 2, 1)
factory.add_circle_arc(2, 3, 4, 2)
factory.add_line(4, 5, 3)
factory.add_circle_arc(5, 6, 7, 4)
factory.add_line(7, 8, 5)
factory.add_line(8, 9, 6)
factory.add_line(9, 10, 7)
factory.add_circle_arc(10, 11, 12, 8)
factory.add_line(12, 13, 9)
factory.add_circle_arc(13, 14, 15, 10)
factory.add_line(15, 16, 11)
factory.add_line(16, 1, 12)

factory.synchronize()

# Физические группы (границы и объем)
factory.add_physical_group(1, [12], -1, "left")
factory.add_physical_group(0, [16, 1], -1, "left")
factory.add_physical_group(1, [6], -1, "right")
factory.add_physical_group(0, [8, 9], -1, "right")


factory.synchronize()
factory.add_curve_loop(collect(1:12), 13)
factory.add_plane_surface([13], 6)
factory.synchronize()
factory.add_physical_group(2, [6], -1, "test_sample")
factory.synchronize()

function add_field_ball!(x, y, r, lc)
  ball_id = gmsh.model.mesh.field.add("Ball")
  gmsh.model.mesh.field.setNumber(ball_id, "XCenter", x)
  gmsh.model.mesh.field.setNumber(ball_id, "YCenter", y)
  gmsh.model.mesh.field.setNumber(ball_id, "ZCenter", 0.0)
  gmsh.model.mesh.field.setNumber(ball_id, "Radius", r)
  gmsh.model.mesh.field.setNumber(ball_id, "Thickness", r * 2)
  gmsh.model.mesh.field.setNumber(ball_id, "VIn", lc)
  gmsh.model.mesh.field.setNumber(ball_id, "VOut", 1e10)
  return ball_id
end

ball_field1 = add_field_ball!(h1 + r, 0, 1.3r, lc/10)
ball_field2 = add_field_ball!(h1 + r + l, 0, 1.3r, lc/10)
ball_field3 = add_field_ball!(h1 + r + l, b_width, 1.3r, lc/10)
ball_field4 = add_field_ball!(h1 + r, b_width, 1.3r, lc/10)

min_id = gmsh.model.mesh.field.add("Min")
gmsh.model.mesh.field.setNumbers(min_id, "FieldsList", [ball_field1, ball_field2, ball_field3, ball_field4])

gmsh.model.mesh.field.setAsBackgroundMesh(min_id)
gmsh.option.setNumber("Mesh.MeshSizeExtendFromBoundary", 0)
gmsh.option.setNumber("Mesh.MeshSizeFromPoints", 0)

gmsh.option.setNumber("Mesh.MeshSizeMax", lc) 

gmsh.model.mesh.generate(2)
model_file = "test_sample.msh" |> tcf
gmsh.write(model_file)
gmsh.finalize()

Геометрия плоского образца со сгенерированной сеткой

2.3.2.2 Матрица жесткости и кинематические связи

Здесь мы задаем физические свойства материала, формируем матрицу упругости для плоского напряженного состояния и определяем функции-операторы (через нотацию Фойгта) для отложенных вычислений в Gridap.

Список 2.2

young_modulus = 2.1e11
poisson_ratio = 0.3
plain_stress = true

k1 = poisson_ratio / (1 - poisson_ratio)
k2 = 0.5 * (1 - 2 * poisson_ratio) / (1 - poisson_ratio)
k3 = (1 - poisson_ratio) / 2

elasticity_matrix = Matrix{Float64}(undef, 3, 3)
if plain_stress
  elasticity_matrix[:, :] = [
    1 k1 0
    k1 1 0
    0 0 k2
  ] * young_modulus * (1 - poisson_ratio) /
    ((1 + poisson_ratio) * (1 - 2 * poisson_ratio))
else
  elasticity_matrix[:, :] = [
    1 poisson_ratio 0
    poisson_ratio 1 0
    0 0 k3
  ] * young_modulus / (1 - poisson_ratio^2)
end

# Преобразование градиента перемещений в вектор деформаций Фойгта
function voigt_strain(grad_u)
  eps_x = grad_u[1, 1]
  eps_y = grad_u[2, 2]
  gamma_xy = grad_u[1, 2] + grad_u[2, 1]
  return VectorValue(eps_x, eps_y, gamma_xy)
end

# Вычисление вектора напряжений по вектору деформаций
function voigt_stress(eps)
  sigma_x = elasticity_matrix[1, 1] * eps[1] +
            elasticity_matrix[1, 2] * eps[2] +
            elasticity_matrix[1, 3] * eps[3]
  sigma_y = elasticity_matrix[2, 1] * eps[1] +
            elasticity_matrix[2, 2] * eps[2] +
            elasticity_matrix[2, 3] * eps[3]
  tau_xy = elasticity_matrix[3, 1] * eps[1] +
           elasticity_matrix[3, 2] * eps[2] +
           elasticity_matrix[3, 3] * eps[3]
  return VectorValue(sigma_x, sigma_y, tau_xy)
end

# Ленивые операторы для интеграции по области
ε(u) = Operation(voigt_strain)(∇(u))
σ(u) = Operation(voigt_stress)(ε(u))

2.3.2.3 Пространства конечных элементов и слабая форма

Здесь определяем векторное двумерное поле и тестовые функции для записи слабой формы, также в модели задаются кинематические граничные условия: на левой границе нулевые смещения, а на правой задается смещение 1мм в положительном направлении оси Х глобальной системы координат (ГСК). Учитывается также, что в модели нет объемных сил.

Список 2.3

# Поле перемещений векторное (2D)
ref_fe = ReferenceFE(lagrangian, VectorValue{2, Float64}, 1)

# Тестовые функции с условиями Дирихле на левом и правом торцах
test_space = TestFESpace(
  model,
  ref_fe,
  conformity = :H1,
  dirichlet_tags = ["left", "right"]
)

# Граничные условия: слева защемление, справа смещение 1 мм
disp_left = VectorValue(0.0, 0.0)
disp_right = VectorValue(1e-3, 0.0)
trial_space = TrialFESpace(test_space, [disp_left, disp_right])

# Мера интегрирования
degree = 2
triangulation = Triangulation(model)
d_omega = Measure(triangulation, degree)

# Слабая форма: билинейная форма (тензорное скалярное произведение)
function weak_form(u, v)
  return ∫( ε(v) ⋅ σ(u) ) * d_omega
end

# Линейная форма равна нулю (объемных сил нет)
function linear_form(v)
  return 0
end

2.3.2.4 Решение задачи

На этом этапе происходит сборка (ансамблирование) глобальной матрицы жесткости и решение СЛАУ метода конечных элементов, используя встроенный решатель BackslashSolver:

Список 2.4

fe_operator = AffineFEOperator(weak_form, linear_form, trial_space, test_space)

initial_guess = zeros(Float64, num_free_dofs(test_space))
solution = FEFunction(trial_space, initial_guess)

linear_solver = BackslashSolver()
fe_solver = LinearFESolver(linear_solver)

# Запуск расчета
solution, _ = solve!(solution, fe_solver, fe_operator)

2.3.3 Визуализация результатов

Для анализа полученных полей деформаций и напряжений мы построим поля перемещений $u_x$ и продольных напряжений $\sigma_x$.

Список 2.5: Визуализация результатов вычислений

fig = Figure(size=(800, 350))

# 1. Поле перемещений
ax1 = Axis(fig[1, 1], aspect=DataAspect(), title="а)", titlealign=:left)
plt1 = warpedmeshtri!(ax1, triangulation, warp=solution, warpscale=20, colormap=:viridis, nomesh=true)
wireframe!(ax1, triangulation, overdraw = true, color=:black, linewidth=0.5, alpha=0.1)

Colorbar(fig[1, 2], plt1, height = Relative(1), label="|u|, m")

# 2. Поле продольных напряжений σ_x
get_sigma_x(s) = s[1]
sigma_x_field = (x) -> (Operation(get_sigma_x)(σ(solution)))(x) / 1e6

ax2 = Axis(fig[2, 1], aspect=DataAspect(), title="б)", titlealign=:left)
plt2 = warpedmeshtri!(ax2, triangulation, warp=solution, warpscale=20, colorby=sigma_x_field, colormap=:coolwarm, nomesh=true)
# wireframe!(ax2, triangulation, overdraw = true, color=:black, linewidth=0.5, alpha=0.1)
Colorbar(fig[2, 2], plt2, height = Relative(1), label="σₓ, МПа")

fig

Рисунок 2.1: Поля деформаций (а) и напряжений (б) образца

2.4 Анализ и интерпретация результатов

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

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

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

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

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

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

Измените кинематические граничные условия: вместо заданного перемещения на правом торце приложите распределенную растягивающую нагрузку (условие Неймана). Сравните полученное поле напряжений.
Выведите график распределения напряжений $\sigma_x$ вдоль центральной продольной оси образца. Оцените коэффициент концентрации напряжений в зоне скруглений.
Исследуйте сходимость решения. Измельчите сетку (параметр lc в Gmsh) и постройте график зависимости максимального значения $\sigma_x$ от характерного размера элемента.
Замените материал образца на алюминиевый сплав (E = 7·10¹⁰ Па, ν = 0.34). Проведите расчет и проанализируйте, как изменились поля перемещений и напряжений по сравнению со стальным образцом при тех же граничных условиях.

```{julia} #| echo: false #| output: false tcf(filename::String) = (@__DIR__) * "\\" * filename; import Pkg; Pkg.activate((@__DIR__) * "\\..\\.."); using Gmsh using Gridap using GridapGmsh using GridapMakie using GLMakie using CCMechBook using Gridap.Arrays using Gridap.Fields using LinearAlgebra ``` # Моделирование растяжения плоского образца ## Аннотация - **Физическая задача**: Моделирование упругого деформирования плоского образца при одноосном растяжении с использованием метода конечных элементов. - **Освоенные численные методы**: Метод Галеркина, конечно-элементная аппроксимация векторных полей, построение матрицы жесткости для плоского напряженного и плоского деформированного состояний. - **Предшествующие знания**: Основы метода конечных элементов, слабая формулировка краевых задач, закон Гука для изотропного материала. ## Основные теоретические сведения ### Упругое деформирование Любые элементы конструкций деформируются под действием приложенных к ним сил. В инженерных расчетах поведение материалов традиционно разделяют на упругое и пластическое. При деформировании внешние силы совершают определенную работу. Если материал ведет себя упруго, эта работа запасается в виде потенциальной энергии деформации, которая полностью высвобождается при снятии нагрузки, возвращая тело в исходное состояние. В случае пластического деформирования часть работы затрачивается на необратимое изменение структуры материала и рассеивается в виде тепла. Чем выше жесткость материала, тем меньше запасенная упругая энергия при одной и той же системе приложенных сил. В пределах малых деформаций зависимость между напряжением и деформацией считается линейной и описывается законом Гука. В обобщенной форме для изотропного материала он записывается следующим образом: $$ \begin{cases} \varepsilon_x = \frac{1}{E}\left[\sigma_x - \nu (\sigma_y + \sigma_z)\right]; \\ \varepsilon_y = \frac{1}{E}\left[\sigma_y - \nu (\sigma_x + \sigma_z)\right]; \\ \varepsilon_z = \frac{1}{E}\left[\sigma_z - \nu (\sigma_x + \sigma_y)\right]; \\ \gamma_{xy} = \frac{\tau_{xy}}{G}; \\ \gamma_{yz} = \frac{\tau_{yz}}{G}; \\ \gamma_{zx} = \frac{\tau_{zx}}{G}; \end{cases} $$ где $\varepsilon$ — линейная деформация вдоль оси; $\sigma$ — нормальное напряжение, Па; $E$ — модуль Юнга, Па; $\nu$ — коэффициент Пуассона, б.р.; $\tau$ — касательное напряжение, Па; $\gamma$ — угловая деформация, б.р.; $G$ — модуль сдвига, $G = \frac{E}{2(1+\nu)}$, Па. **Коэффициент Пуассона** показывает, во сколько раз изменяется поперечное сечение деформированного тела при его растяжении или сжатии (отношение относительного поперечного сжатия к относительному продольному растяжению) и зависит от природы материала. **Модуль Юнга** (модуль упругости) характеризует сопротивление материала растяжению и сжатию при упругой деформации. Эти два коэффициента полностью описывают упругие свойства изотропного материала. ### Нотация Фойгта Закон Гука выше можно представить в виде линейного отображения: $$ \sigma = {}^4\mathbf{C}_{\mathrm{eff}}{:}\varepsilon $$ где $\sigma$ и $\varepsilon$ - тензоры напряжений и деформаций, а ${}^4\mathbf{C}$ - тензор упругих постоянных. Поскольку напряжения и деформации являются тензорами второго порядка, отсюда следует что ${}^4\mathbf{C}$ является тензором четвертого порядка. Однако в силу того, что такое представление не очень удобно для записи, широкое распространение получила нотация Фойгта. Согласно этой нотации, если учесть симметрию тензоров напряжений и деформаций, то их можно представить в виде векторов, содержащих 6 компонент: $$ \sigma_v = \begin{Bmatrix} \sigma_{11} \\ \sigma_{22} \\ \sigma_{33} \\ \sigma_{23} \\ \sigma_{13} \\ \sigma_{12} \end{Bmatrix}= \begin{Bmatrix} \sigma_x \\ \sigma_y \\ \sigma_z \\ \tau_{yz} \\ \tau_{xz} \\ \tau_{xy} \end{Bmatrix} \qquad \varepsilon_v = \begin{Bmatrix} \varepsilon_{11} \\ \varepsilon_{22} \\ \varepsilon_{33} \\ 2\varepsilon_{23} \\ 2\varepsilon_{13} \\ 2\varepsilon_{12} \end{Bmatrix}= \begin{Bmatrix} \varepsilon_x \\ \varepsilon_y \\ \varepsilon_z \\ \gamma_{yz} \\ \gamma_{xz} \\ \gamma_{xy} \end{Bmatrix} $$ Тогда тензор упругих постоянных примет вид матрицы: $$ \mathbf{C}_v = \frac{E(1-\nu)}{(1+\nu)(1-2\nu)} \begin{bmatrix} 1 & \frac{\nu}{1-\nu} & \frac{\nu}{1-\nu} & 0 & 0 & 0 \\ \frac{\nu}{1-\nu} & 1 & \frac{\nu}{1-\nu} & 0 & 0 & 0 \\ \frac{\nu}{1-\nu} & \frac{\nu}{1-\nu} & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{1-2\nu}{2(1-\nu)} & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{1-2\nu}{2(1-\nu)} & 0 \\ 0 & 0 & 0 & 0 & 0 & \frac{1-2\nu}{2(1-\nu)} \end{bmatrix} $$ Сам закон Гука в матричном виде записывается компактно: $$ \sigma_v = \mathbf{C}_v \varepsilon_v $$ ### Плоское напряженное и плоское деформированное состояния Приведенную выше зависимость можно существенно упростить, если рассматриваемая область является плоской, а нагрузки или смещения действуют только в плоскости контура. В этом случае выделяют два классических варианта: 1. **Плоское напряженное состояние** — напряжения вдоль третьей оси полагаются нулевыми, а деформации свободны (тело имеет малую толщину, например, тонкая пластина). 2. **Плоское деформированное состояние** — деформации (и перемещения) вдоль третьей оси жестко зафиксированы, а напряжения остаются свободными (тело имеет большую толщину, например, протяженная плотина). Для плоского деформированного состояния матрица упругости имеет вид: $$ \mathbf{C}_v = \frac{E(1-\nu)}{(1+\nu)(1-2\nu)} \begin{bmatrix} 1 & \frac{\nu}{1-\nu} & 0 \\ \frac{\nu}{1-\nu} & 1 & 0 \\ 0 & 0 & \frac{1-2\nu}{2(1-\nu)} \end{bmatrix} $$ Для плоского напряженного состояния: $$ \mathbf{C}_{v} = \frac{E}{1-\nu^2} \begin{bmatrix} 1 & \nu & 0 \\ \nu & 1 & 0 \\ 0 & 0 & \frac{1-\nu}{2} \end{bmatrix} $$ ### Слабая формулировка задачи упругости Рассмотрим упругое тело, занимающее область $\Omega$. В состоянии статического равновесия для него справедливо уравнение Навье: $$ -\nabla \cdot \sigma(u) = f \quad \text{в} \quad \Omega $$ где $\sigma(u)$ — тензор напряжений, $f$ — вектор объемных сил. Для получения слабой формы умножим обе части уравнения на произвольную векторную пробную функцию $v$ и проинтегрируем по объему $\Omega$: $$ -\int_{\Omega} (\nabla \cdot \sigma(u)) \cdot v \, dx = \int_{\Omega} f \cdot v \, dx $$ Применяя интегрирование по частям (формулу Гаусса-Остроградского) к члену с дивергенцией, получаем: $$ -\int_{\Omega} (\nabla \cdot \sigma(u)) \cdot v \, dx = \int_{\Omega} \sigma(u) \cdot \nabla v \, dx - \int_{\partial\Omega} (\sigma(u) \cdot n) \cdot v \, dx $$ где $n$ — внешняя нормаль к границе $\partial\Omega$. Величина $\sigma(u) \cdot n$ представляет собой вектор поверхностных сил (напряжений) на границе. Так как тензор напряжений $\sigma$ симметричен, его внутреннее произведение с антисимметричной частью градиента пробной функции $\nabla v$ равно нулю. Поэтому полный градиент $\nabla v$ можно заменить на его симметричную часть $\varepsilon(v)$ (тензор деформаций пробного поля). Это приводит к билинейной форме: $$ a(u, v) = \int_{\Omega} \sigma(u) \cdot \varepsilon(v) \, dx $$ Линейная форма (правая часть уравнения) с учетом граничных условий Неймана $\Gamma_N$ примет вид: $$ L(v) = \int_{\Omega} f \cdot v \, dx + \int_{\Gamma_N} \bar{t} \cdot v \, ds $$ где $\bar{t}$ — заданный вектор поверхностных нагрузок. В нашем случае мы будем рассматривать задачу без объемных сил ($f = 0$) и без распределенных поверхностных нагрузок ($\bar{t} = 0$), применяя только кинематические граничные условия (перемещения). ## Практическая задача ### Постановка задачи Смоделируем растяжение плоского стального образца. Геометрия соответствует стандартному плоскому образцу по ГОСТ 1497-84 (тип I, №6). На левой границе образца перемещения жестко зафиксированы (защемление), а правая граница смещается на **1 мм** вдоль продольной оси. **Исходные данные:** * Материал: Сталь * Модуль Юнга: **E = 2.1·10¹¹ Па** * Коэффициент Пуассона: **ν = 0,3** * Толщина образца: **B = 0,04 м** * Начальная длина: **l₀ = 0,14 м** ### Реализация в программном стеке #### Создание геометрии в Gmsh В этом блоке мы формируем параметрическую 2D модель плоского образца со скруглениями (галтелями) и генерируем сетку. ```{julia} #| label: lst-geometry-creation #| lst-cap: "Создание геометрии и сетки плоского образца в Gmsh" #| output: false order = 1 gmsh.initialize() b_width = 0.04 h1 = 0.08 r = (b_width - 0.03) / 2 l = 0.14 + 2 * sqrt(0.02 * 0.03) lc = 0.005 gmsh.option.set_number("General.Terminal", 1) gmsh.option.set_number("Mesh.ElementOrder", order) gmsh.option.set_number("Mesh.Algorithm", 5) gmsh.model.add("flat_sample") factory = gmsh.model.geo # Точки контура factory.add_point(0, 0, 0, lc, 1) factory.add_point(h1, 0, 0, lc, 2) factory.add_point(h1 + r, 0, 0, lc, 3) factory.add_point(h1 + r, r, 0, lc, 4) factory.add_point(h1 + r + l, r, 0, lc, 5) factory.add_point(h1 + r + l, 0, 0, lc, 6) factory.add_point(h1 + r + l + r, 0, 0, lc, 7) factory.add_point(h1 + r + l + r + h1, 0, 0, lc, 8) factory.add_point(h1 + r + l + r + h1, b_width, 0, lc, 9) factory.add_point(h1 + r + l + r, b_width, 0, lc, 10) factory.add_point(h1 + r + l, b_width, 0, lc, 11) factory.add_point(h1 + r + l, b_width - r, 0, lc, 12) factory.add_point(h1 + r, b_width - r, 0, lc, 13) factory.add_point(h1 + r, b_width, 0, lc, 14) factory.add_point(h1, b_width, 0, lc, 15) factory.add_point(0, b_width, 0, lc, 16) # Линии и дуги factory.add_line(1, 2, 1) factory.add_circle_arc(2, 3, 4, 2) factory.add_line(4, 5, 3) factory.add_circle_arc(5, 6, 7, 4) factory.add_line(7, 8, 5) factory.add_line(8, 9, 6) factory.add_line(9, 10, 7) factory.add_circle_arc(10, 11, 12, 8) factory.add_line(12, 13, 9) factory.add_circle_arc(13, 14, 15, 10) factory.add_line(15, 16, 11) factory.add_line(16, 1, 12) factory.synchronize() # Физические группы (границы и объем) factory.add_physical_group(1, [12], -1, "left") factory.add_physical_group(0, [16, 1], -1, "left") factory.add_physical_group(1, [6], -1, "right") factory.add_physical_group(0, [8, 9], -1, "right") factory.synchronize() factory.add_curve_loop(collect(1:12), 13) factory.add_plane_surface([13], 6) factory.synchronize() factory.add_physical_group(2, [6], -1, "test_sample") factory.synchronize() function add_field_ball!(x, y, r, lc) ball_id = gmsh.model.mesh.field.add("Ball") gmsh.model.mesh.field.setNumber(ball_id, "XCenter", x) gmsh.model.mesh.field.setNumber(ball_id, "YCenter", y) gmsh.model.mesh.field.setNumber(ball_id, "ZCenter", 0.0) gmsh.model.mesh.field.setNumber(ball_id, "Radius", r) gmsh.model.mesh.field.setNumber(ball_id, "Thickness", r * 2) gmsh.model.mesh.field.setNumber(ball_id, "VIn", lc) gmsh.model.mesh.field.setNumber(ball_id, "VOut", 1e10) return ball_id end ball_field1 = add_field_ball!(h1 + r, 0, 1.3r, lc/10) ball_field2 = add_field_ball!(h1 + r + l, 0, 1.3r, lc/10) ball_field3 = add_field_ball!(h1 + r + l, b_width, 1.3r, lc/10) ball_field4 = add_field_ball!(h1 + r, b_width, 1.3r, lc/10) min_id = gmsh.model.mesh.field.add("Min") gmsh.model.mesh.field.setNumbers(min_id, "FieldsList", [ball_field1, ball_field2, ball_field3, ball_field4]) gmsh.model.mesh.field.setAsBackgroundMesh(min_id) gmsh.option.setNumber("Mesh.MeshSizeExtendFromBoundary", 0) gmsh.option.setNumber("Mesh.MeshSizeFromPoints", 0) gmsh.option.setNumber("Mesh.MeshSizeMax", lc) gmsh.model.mesh.generate(2) model_file = "test_sample.msh" |> tcf gmsh.write(model_file) gmsh.finalize() ``` ```{julia} #| label: mesh #| fig-cap: "Геометрия плоского образца со сгенерированной сеткой" #| output: true #| echo: false model = GmshDiscreteModelNoVerbose(model_file); Ω = Triangulation(model) Γn = BoundaryTriangulation(model, tags=["left"]); Γd = BoundaryTriangulation(model, tags=["right"]); fig = Figure(size=(800, 160)); ax = Axis(fig[1,1], aspect=DataAspect()); plt = plot!(ax, Ω); wireframe!(Ω, color=:black, linewidth=0.5); plot!(Γn, color=:blue, linewidth=2); plot!(Γd, color=:red, linewidth=2); fig ``` #### Матрица жесткости и кинематические связи Здесь мы задаем физические свойства материала, формируем матрицу упругости для плоского напряженного состояния и определяем функции-операторы (через нотацию Фойгта) для отложенных вычислений в Gridap. ```{julia} #| label: lst-elasticity-constants #| lst-cap: "Определение констант упругости и операторов Фойгта" #| output: false young_modulus = 2.1e11 poisson_ratio = 0.3 plain_stress = true k1 = poisson_ratio / (1 - poisson_ratio) k2 = 0.5 * (1 - 2 * poisson_ratio) / (1 - poisson_ratio) k3 = (1 - poisson_ratio) / 2 elasticity_matrix = Matrix{Float64}(undef, 3, 3) if plain_stress elasticity_matrix[:, :] = [ 1 k1 0 k1 1 0 0 0 k2 ] * young_modulus * (1 - poisson_ratio) / ((1 + poisson_ratio) * (1 - 2 * poisson_ratio)) else elasticity_matrix[:, :] = [ 1 poisson_ratio 0 poisson_ratio 1 0 0 0 k3 ] * young_modulus / (1 - poisson_ratio^2) end # Преобразование градиента перемещений в вектор деформаций Фойгта function voigt_strain(grad_u) eps_x = grad_u[1, 1] eps_y = grad_u[2, 2] gamma_xy = grad_u[1, 2] + grad_u[2, 1] return VectorValue(eps_x, eps_y, gamma_xy) end # Вычисление вектора напряжений по вектору деформаций function voigt_stress(eps) sigma_x = elasticity_matrix[1, 1] * eps[1] + elasticity_matrix[1, 2] * eps[2] + elasticity_matrix[1, 3] * eps[3] sigma_y = elasticity_matrix[2, 1] * eps[1] + elasticity_matrix[2, 2] * eps[2] + elasticity_matrix[2, 3] * eps[3] tau_xy = elasticity_matrix[3, 1] * eps[1] + elasticity_matrix[3, 2] * eps[2] + elasticity_matrix[3, 3] * eps[3] return VectorValue(sigma_x, sigma_y, tau_xy) end # Ленивые операторы для интеграции по области ε(u) = Operation(voigt_strain)(∇(u)) σ(u) = Operation(voigt_stress)(ε(u)) ``` #### Пространства конечных элементов и слабая форма Здесь определяем векторное двумерное поле и тестовые функции для записи слабой формы, также в модели задаются кинематические граничные условия: на левой границе нулевые смещения, а на правой задается смещение 1мм в положительном направлении оси Х глобальной системы координат (ГСК). Учитывается также, что в модели нет объемных сил. ```{julia} #| label: lst-weak-form-elasticity #| lst-cap: "Инициализация пространств КЭ и слабой формы" #| output: false # Поле перемещений векторное (2D) ref_fe = ReferenceFE(lagrangian, VectorValue{2, Float64}, 1) # Тестовые функции с условиями Дирихле на левом и правом торцах test_space = TestFESpace( model, ref_fe, conformity = :H1, dirichlet_tags = ["left", "right"] ) # Граничные условия: слева защемление, справа смещение 1 мм disp_left = VectorValue(0.0, 0.0) disp_right = VectorValue(1e-3, 0.0) trial_space = TrialFESpace(test_space, [disp_left, disp_right]) # Мера интегрирования degree = 2 triangulation = Triangulation(model) d_omega = Measure(triangulation, degree) # Слабая форма: билинейная форма (тензорное скалярное произведение) function weak_form(u, v) return ∫( ε(v) ⋅ σ(u) ) * d_omega end # Линейная форма равна нулю (объемных сил нет) function linear_form(v) return 0 end ``` #### Решение задачи На этом этапе происходит сборка (ансамблирование) глобальной матрицы жесткости и решение СЛАУ метода конечных элементов, используя встроенный решатель `BackslashSolver`: ```{julia} #| label: lst-solve-elasticity #| lst-cap: "Сборка СЛАУ и запуск решателя" #| output: false fe_operator = AffineFEOperator(weak_form, linear_form, trial_space, test_space) initial_guess = zeros(Float64, num_free_dofs(test_space)) solution = FEFunction(trial_space, initial_guess) linear_solver = BackslashSolver() fe_solver = LinearFESolver(linear_solver) # Запуск расчета solution, _ = solve!(solution, fe_solver, fe_operator) ``` ### Визуализация результатов Для анализа полученных полей деформаций и напряжений мы построим поля перемещений $u_x$ и продольных напряжений $\sigma_x$. ```{julia} #| label: fig-vis-elasticity #| lst-label: lst-vis-elasticity #| lst-cap: "Визуализация результатов вычислений" #| fig-cap: "Поля деформаций (*а*) и напряжений (*б*) образца" #| output: true #| fig = Figure(size=(800, 350)) # 1. Поле перемещений ax1 = Axis(fig[1, 1], aspect=DataAspect(), title="а)", titlealign=:left) plt1 = warpedmeshtri!(ax1, triangulation, warp=solution, warpscale=20, colormap=:viridis, nomesh=true) wireframe!(ax1, triangulation, overdraw = true, color=:black, linewidth=0.5, alpha=0.1) Colorbar(fig[1, 2], plt1, height = Relative(1), label="|u|, m") # 2. Поле продольных напряжений σ_x get_sigma_x(s) = s[1] sigma_x_field = (x) -> (Operation(get_sigma_x)(σ(solution)))(x) / 1e6 ax2 = Axis(fig[2, 1], aspect=DataAspect(), title="б)", titlealign=:left) plt2 = warpedmeshtri!(ax2, triangulation, warp=solution, warpscale=20, colorby=sigma_x_field, colormap=:coolwarm, nomesh=true) # wireframe!(ax2, triangulation, overdraw = true, color=:black, linewidth=0.5, alpha=0.1) Colorbar(fig[2, 2], plt2, height = Relative(1), label="σₓ, МПа") fig ``` ## Анализ и интерпретация результатов Анализ полученных полей показывает, что максимальные напряжения концентрируются в зонах галтелей (переходов от широкой части к узкой). Это классический пример концентрации напряжений, что полностью соответствует аналитической теории упругости. В центральной (рабочей) части образца напряжения распределены равномерно, обеспечивая корректные условия чистого одноосного растяжения. Перемещения вдоль продольной оси линейно возрастают от левой закрепленной границы (0 мм) до правой смещаемой (1 мм). Использование масштабного коэффициента при визуализации деформированного состояния (график 3) позволяет инженеру наглядно оценить характер искажения геометрии, включая поперечное сужение (эффект Пуассона). ## Контрольные вопросы и задания ### Теоретические вопросы 1. Какой физический смысл имеет коэффициент Пуассона и каковы его предельные значения для реальных изотропных материалов? 2. В чем заключается фундаментальное различие между плоским напряженным и плоским деформированным состояниями? Приведите примеры конструкций для каждого случая. 3. Почему в билинейной (слабой) форме задачи упругости используется симметричная часть градиента (тензор деформаций), а не полный градиент перемещений? 4. Какие граничные условия необходимы и достаточны для исключения перемещения тела как жесткого целого? ### Практические задания 1. Измените кинематические граничные условия: вместо заданного перемещения на правом торце приложите распределенную растягивающую нагрузку (условие Неймана). Сравните полученное поле напряжений. 2. Выведите график распределения напряжений $\sigma_x$ вдоль центральной продольной оси образца. Оцените коэффициент концентрации напряжений в зоне скруглений. 3. Исследуйте сходимость решения. Измельчите сетку (параметр `lc` в `Gmsh`) и постройте график зависимости максимального значения $\sigma_x$ от характерного размера элемента. 4. Замените материал образца на алюминиевый сплав (**E = 7·10¹⁰ Па**, **ν = 0.34**). Проведите расчет и проанализируйте, как изменились поля перемещений и напряжений по сравнению со стальным образцом при тех же граничных условиях.