3 Расчет напряженно-деформированного состояния гибкого шарнира при больших перемещениях и поворотах

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

Предварительные знания: Основы линейной теории упругости (Урок 2), тензорное исчисление, базовые навыки работы с Gmsh.
Освоенные численные методы: Импорт геометрии из форматов САПР (STEP), формулировка задачи с учетом больших деформаций (тензор деформаций Грина-Лагранжа, 2-й тензор напряжений Пиолы-Кирхгофа), метод Ньютона-Рафсона для нелинейной МКЭ-задачи, инкрементальное нагружение.
Физическая задача: Расчет напряженно-деформированного состояния (НДС) трехмерного шарнирного узла при больших перемещениях и поворотах (геометрическая нелинейность).

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

3.2.1 Геометрически нелинейная упругость

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

Для описания больших деформаций вводится градиент деформации $\mathbf{F}$, связывающий координаты точек в начальной и деформированной конфигурациях:

\[ \mathbf{F} = \mathbf{I} + \nabla \mathbf{u} \tag{3.1}\]

где $\mathbf{I}$ — единичный тензор, $\nabla \mathbf{u}$ — тензор градиента перемещений.

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

\[ \mathbf{E} = \frac{1}{2}(\mathbf{F}^T \mathbf{F} - \mathbf{I}) \tag{3.2}\]

Связь между напряжениями и деформациями в нелинейной постановке требует использования объективных мер напряжений. Одной из таких мер является второй тензор напряжений Пиолы-Кирхгофа $\mathbf{S}$. Для изотропного гиперупругого материала, подчиняющегося модели Сен-Венана — Кирхгофа (расширение закона Гука на большие деформации), определяющее соотношение имеет вид:

\[ \mathbf{S}(\mathbf{u}) = \lambda \mathrm{tr}(\mathbf{E})\mathbf{I} + 2\mu\mathbf{E} \tag{3.3}\]

где $\lambda$ и $\mu$ — параметры Ламе. Модель Сен-Венана — Кирхгофа хорошо работает для больших перемещений и поворотов, но применима только при малых относительных удлинениях.

3.2.2 Нелинейная слабая форма и линеаризация

Поскольку тензор деформаций $\mathbf{E}$ нелинейно зависит от перемещений $\mathbf{u}$, уравнение равновесия становится нелинейным относительно $\mathbf{u}$. Построение слабой формы (уравнения баланса виртуальных работ) дает невязку (residual) $\mathcal{R}(\mathbf{u}, \mathbf{v})$, которая в состоянии равновесия должна равняться нулю:

\[ \mathcal{R}(\mathbf{u}, \mathbf{v}) = \int_\Omega \delta\mathbf{E}(\mathbf{v}, \mathbf{u}) : \mathbf{S}(\mathbf{u}) \, d\Omega - \int_\Gamma \mathbf{v} \cdot \mathbf{t} \, d\Gamma = 0 \tag{3.4}\]

где $\mathbf{v}$ — векторная пробная функция (виртуальные перемещения), $\mathbf{t}$ — поверхностная нагрузка, а вариация деформаций Грина-Лагранжа определяется как:

\[ \delta\mathbf{E}(\mathbf{v}, \mathbf{u}) = \frac{1}{2}\left( \nabla \mathbf{v} \cdot \mathbf{F}(\mathbf{u}) + \mathbf{F}(\mathbf{u})^T \cdot \nabla \mathbf{v}^T \right) \tag{3.5}\]

Для решения нелинейного уравнения $\mathcal{R}(\mathbf{u}, \mathbf{v}) = 0$ применяется итерационный метод Ньютона-Рафсона. Он требует вычисления производной Фреше по направлению $\Delta \mathbf{u}$ — Якобиана (матрицы жесткости). В геометрически нелинейных задачах Якобиан состоит из двух слагаемых: материальной и геометрической матриц жесткости:

\[ \mathcal{J}(\mathbf{u}, \Delta\mathbf{u}, \mathbf{v}) = \int_\Omega \delta\mathbf{E}(\mathbf{v}, \mathbf{u}) : d\mathbf{S}(\Delta\mathbf{u}, \mathbf{u}) \, d\Omega + \int_\Omega \nabla \mathbf{v} : (\nabla \Delta\mathbf{u} \cdot \mathbf{S}(\mathbf{u})) \, d\Omega \tag{3.6}\]

Итерационный процесс продолжается до тех пор, пока невязка не станет меньше заданного допуска.

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

3.3.1 Постановка задачи и импорт геометрии

Мы рассмотрим стальной кронштейн (шарнир), левый торец которого жестко заделан, а на правый наложено кинематическое граничное условие — поворот на заданный угол и смещение. Вместо ручного построения геометрии мы импортируем готовую модель из файла САПР model.step.

Список 3.1: Импорт CAD геометрии и построение сетки

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

# Импорт геометрии из STEP файла
path = "model.step" |> tcf
gmsh.model.occ.importShapes(path)
gmsh.model.occ.synchronize()

# Получение 3D объектов и их габаритов для автоматической разметки
d3, = gmsh.model.getEntities(3)
xmin, ymin, zmin, xmax, ymax, zmax = gmsh.model.occ.getBoundingBox(d3...)
EPS = 1e-5

# Поиск поверхностей, лежащих на левом и правом краях по координате Z
left_bnd = gmsh.model.getEntitiesInBoundingBox(-Inf, -Inf, zmin - EPS, +Inf, +Inf, zmin + EPS)
right_bnd = gmsh.model.getEntitiesInBoundingBox(-Inf, -Inf, zmax - EPS, +Inf, +Inf, zmax + EPS)

# Назначение физических групп (0, 1, 2 - точки, линии, поверхности)
gmsh.model.addPhysicalGroup(3, [d3[2]], -1, "domain")
for dim in (0,1,2)
    gmsh.model.addPhysicalGroup(dim, [tag[2] for tag in left_bnd if tag[1] == dim], -1, "fixed")
    gmsh.model.addPhysicalGroup(dim, [tag[2] for tag in right_bnd if tag[1] == dim], -1, "loaded")
end

gmsh.model.occ.synchronize()

# Настройки генератора сетки
gmsh.option.setNumber("Mesh.Algorithm", 4)
gmsh.option.setNumber("Mesh.MeshSizeMax", 5.0)

gmsh.model.mesh.generate(3)
msh_file = "hinge_model.msh" |> tcf
gmsh.write(msh_file)
gmsh.finalize()

3.3.2 Дискретизация и математическая модель в Gridap

Зададим свойства материала (сталь, $E = 2.1\cdot 10^5$ МПа, $\nu = 0.3$) и переведем формулы тензорной кинематики в синтаксис Gridap. Мы будем активно использовать ленивые операторы композиции ∘ и внутреннего произведения ⊙.

Список 3.2: Определение нелинейной кинематики и функционалов

model = GmshDiscreteModel(msh_file)

E_mod = 2.1e5  # Масштабированный модуль
ν = 0.3
λ = (E_mod * ν) / ((1 + ν) * (1 - 2 * ν))
μ = E_mod / (2 * (1 + ν))

# Симметричная часть тензора
sym(v) = 0.5 * (v + v')

# Градиент деформации F
F(∇u) = one(∇u) + ∇u'

# Тензор деформаций Грина-Лагранжа E
Ε(∇u) = 0.5 * (F(∇u)' ⋅ F(∇u) - one(F(∇u)))

# Вариация тензора деформаций dE
dΕ(∇du, ∇u) = sym(∇du ⋅ F(∇u))

# Второй тензор напряжений Пиолы-Кирхгофа S (модель Сен-Венана - Кирхгофа)
S(∇u) = λ * tr(Ε(∇u)) * one(Ε(∇u)) + 2 * μ * Ε(∇u)

# Вариация тензора напряжений dS
dS(∇du, ∇u) = λ * tr(dΕ(∇du, ∇u)) * one(dΕ(∇du, ∇u)) + 2 * μ * dΕ(∇du, ∇u)

degree = 2
Ω = Triangulation(model)
dΩ = Measure(Ω, degree)

# Невязка (residual)
res(u, v) = ∫( (dΕ ∘ (∇(v), ∇(u))) ⊙ (S ∘ ∇(u)) ) * dΩ

# Якобиан: материальная и геометрическая части
jac_mat(u, du, v) = ∫( (dΕ ∘ (∇(v), ∇(u))) ⊙ (dS ∘ (∇(du), ∇(u))) ) * dΩ
jac_geo(u, du, v) = ∫( ∇(v) ⊙ ( (S ∘ ∇(u)) ⋅ ∇(du) ) ) * dΩ
jac(u, du, v) = jac_mat(u, du, v) + jac_geo(u, du, v)

3.3.3 Настройка решателя и инкрементальное нагружение

В нелинейных задачах полное приложение нагрузки за один шаг часто приводит к расходимости метода Ньютона-Рафсона. Чтобы обойти эту проблему, мы разобьем кинематическое нагружение на несколько инкрементальных шагов.

Список 3.3: Инкрементальное решение задачи с использованием метода Ньютона

using LineSearches: BackTracking

reffe = ReferenceFE(lagrangian, VectorValue{3,Float64}, 1)
V = TestFESpace(model, reffe, conformity=:H1, dirichlet_tags=["fixed", "loaded"])

# Настройка нелинейного решателя
nls = NLSolver(show_trace=false, method=:newton, linesearch=BackTracking(), iterations=20)
solver = FESolver(nls)

function solve_step(x0, angle, cache)
    # Нулевое перемещение на "fixed"
    g0 = VectorValue(0.0, 0.0, 0.0)
    
    # Заданное кинематическое вращение и смещение на "loaded"
    g1(x) = VectorValue(
        x[1] * cos(angle) + x[3] * sin(angle) - x[1],
        0.0, 
        -x[1] * sin(angle) + x[3] * cos(angle) - x[3]
    )
    
    U = TrialFESpace(V, [g0, g1])
    op = FEOperator(res, jac, U, V)
    
    uh_step = FEFunction(U, x0)
    uh_step, cache = solve!(uh_step, solver, op, cache)
    
    return uh_step, get_free_dof_values(uh_step), cache
end

disp_max = deg2rad(45) # Поворот на 45 градусов
nsteps = 5
x_vals = zeros(Float64, num_free_dofs(V))
cache = nothing
uh_final = nothing

for step in 1:nsteps
    current_angle = step * disp_max / nsteps
    uh_final, x_vals, cache = solve_step(x_vals, current_angle, cache)
end

3.3.4 Визуализация результатов с помощью Makie

Визуализируем норму вектора перемещений на расчетной модели.

Список 3.4: Визуализация поля перемещений

fig = Figure(size = (600, 600))
ax = Axis3(fig[1,1], aspect = :data)

# Вычисляем норму перемещений
u_mag = sqrt ∘ (uh_final ⋅ uh_final)
u_max = maximum(u_mag.(Ω.grid.node_coordinates));
∂Ω = BoundaryTriangulation(model)
# Отрисовка геометрии с наложением цветовой карты
plt = plot!(ax, ∂Ω, u_mag, colormap=cgrad(:turbo, 10, categorical=true), colorrange=(0,0.99*u_max))
wireframe!(ax, ∂Ω, color=:black, linewidth=1.2)

Colorbar(fig[1, 2], plt, label = "мм")
writevtk(Ω, tcf("disps"); cellfields=Dict("uh"=>uh_final, "u_mag"=>u_mag))
# save("fig.png", fig)
fig

Рисунок 3.1: Распределение нормы перемещений в кронштейне при больших деформациях

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

Использование геометрически нелинейной постановки позволило корректно рассчитать поведение кронштейна при больших поворотах (до 45 градусов). Ключевые выводы:

Линейная модель (малые деформации $\boldsymbol{\varepsilon} = \frac{1}{2}(\nabla \mathbf{u} + \nabla \mathbf{u}^T)$) в такой задаче привела бы к возникновению огромных фиктивных напряжений и нефизичному увеличению объема детали в зоне вращения.

Использование тензора Грина-Лагранжа $\mathbf{E}$ отфильтровывает повороты как жесткого целого, позволяя материалу деформироваться корректно.

Инкрементальный подход (разбиение нагрузки на 5 шагов) — критически важная техника. Матрица Якоби зависит от текущих перемещений $\mathbf{u}$, и при попытке приложить полный угол поворота за одну итерацию алгоритм Ньютона-Рафсона ушел бы в бесконечность из-за плохой начальной догадки.

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

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

В чем физический смысл разницы между градиентом перемещений $\nabla \mathbf{u}$ и тензором Грина-Лагранжа $\mathbf{E}$?

Почему модель материала Сен-Венана — Кирхгофа нельзя применять для эластомеров (резин) при сильном растяжении, несмотря на то, что она допускает большие перемещения?

Каков физический и математический смысл геометрической матрицы жесткости в уравнении Якобиана?

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

Измените количество шагов инкрементального нагружения nsteps до 1. Изучите вывод решателя. Произойдет ли сходимость?
Измените угол поворота disp_max на 60 и 90 градусов. Как изменяется требуемое количество итераций на шаг?
(*) Реализуйте и визуализируйте расчет эквивалентных напряжений по Мизесу (через тензор Коши, который получается путем пуш-форварда тензора Пиолы-Кирхгофа: $\boldsymbol{\sigma} = J^{-1} \mathbf{F} \mathbf{S} \mathbf{F}^T$).

```{julia} #| echo: false #| output: false tcf(filename::String) = (@__DIR__) * "\\" * filename; import Pkg; Pkg.activate((@__DIR__) * "\\..\\.."); using Gmsh using Gridap using GridapGmsh using GLMakie using GridapMakie ``` # Расчет напряженно-деформированного состояния гибкого шарнира при больших перемещениях и поворотах ## Аннотация и цели обучения - **Предварительные знания**: Основы линейной теории упругости (Урок 2), тензорное исчисление, базовые навыки работы с Gmsh. - **Освоенные численные методы**: Импорт геометрии из форматов САПР (STEP), формулировка задачи с учетом больших деформаций (тензор деформаций Грина-Лагранжа, 2-й тензор напряжений Пиолы-Кирхгофа), метод Ньютона-Рафсона для нелинейной МКЭ-задачи, инкрементальное нагружение. - **Физическая задача**: Расчет напряженно-деформированного состояния (НДС) трехмерного шарнирного узла при больших перемещениях и поворотах (геометрическая нелинейность). ## Основные теоретические сведения ### Геометрически нелинейная упругость В предыдущем уроке мы рассматривали классическую линейную теорию упругости, где предполагалось, что перемещения тел бесконечно малы. Это позволяло нам не различать начальную и деформированную конфигурации тела. Однако в реальных инженерных задачах конструкции часто подвергаются значительным изгибам и поворотам, при которых линейная модель дает существенную ошибку (например, искусственно увеличивает объем тела при вращении). Для описания больших деформаций вводится градиент деформации $\mathbf{F}$, связывающий координаты точек в начальной и деформированной конфигурациях: $$ \mathbf{F} = \mathbf{I} + \nabla \mathbf{u} $$ {#eq-deformation-gradient} где $\mathbf{I}$ — единичный тензор, $\nabla \mathbf{u}$ — тензор градиента перемещений. Поскольку градиент деформации содержит в себе поворот тела как жесткого целого (который не должен вызывать напряжений), для оценки истинного изменения формы используется **тензор деформаций Грина-Лагранжа**: $$ \mathbf{E} = \frac{1}{2}(\mathbf{F}^T \mathbf{F} - \mathbf{I}) $$ {#eq-green-lagrange-strain-tensor} Связь между напряжениями и деформациями в нелинейной постановке требует использования объективных мер напряжений. Одной из таких мер является второй **тензор напряжений Пиолы-Кирхгофа** $\mathbf{S}$. Для изотропного гиперупругого материала, подчиняющегося модели Сен-Венана — Кирхгофа (расширение закона Гука на большие деформации), определяющее соотношение имеет вид: $$ \mathbf{S}(\mathbf{u}) = \lambda \mathrm{tr}(\mathbf{E})\mathbf{I} + 2\mu\mathbf{E} $$ {#eq-piola-kirchhoff-stress-tensor} где $\lambda$ и $\mu$ — параметры Ламе. Модель Сен-Венана — Кирхгофа хорошо работает для больших перемещений и поворотов, но применима только при малых относительных удлинениях. ### Нелинейная слабая форма и линеаризация Поскольку тензор деформаций $\mathbf{E}$ нелинейно зависит от перемещений $\mathbf{u}$, уравнение равновесия становится нелинейным относительно $\mathbf{u}$. Построение слабой формы (уравнения баланса виртуальных работ) дает невязку (residual) $\mathcal{R}(\mathbf{u}, \mathbf{v})$, которая в состоянии равновесия должна равняться нулю: $$ \mathcal{R}(\mathbf{u}, \mathbf{v}) = \int_\Omega \delta\mathbf{E}(\mathbf{v}, \mathbf{u}) : \mathbf{S}(\mathbf{u}) \, d\Omega - \int_\Gamma \mathbf{v} \cdot \mathbf{t} \, d\Gamma = 0 $$ {#eq-residual} где $\mathbf{v}$ — векторная пробная функция (виртуальные перемещения), $\mathbf{t}$ — поверхностная нагрузка, а вариация деформаций Грина-Лагранжа определяется как: $$ \delta\mathbf{E}(\mathbf{v}, \mathbf{u}) = \frac{1}{2}\left( \nabla \mathbf{v} \cdot \mathbf{F}(\mathbf{u}) + \mathbf{F}(\mathbf{u})^T \cdot \nabla \mathbf{v}^T \right) $$ {#eq-strain-variation} Для решения нелинейного уравнения $\mathcal{R}(\mathbf{u}, \mathbf{v}) = 0$ применяется итерационный метод Ньютона-Рафсона. Он требует вычисления производной Фреше по направлению $\Delta \mathbf{u}$ — Якобиана (матрицы жесткости). В геометрически нелинейных задачах Якобиан состоит из двух слагаемых: **материальной** и **геометрической** матриц жесткости: $$ \mathcal{J}(\mathbf{u}, \Delta\mathbf{u}, \mathbf{v}) = \int_\Omega \delta\mathbf{E}(\mathbf{v}, \mathbf{u}) : d\mathbf{S}(\Delta\mathbf{u}, \mathbf{u}) \, d\Omega + \int_\Omega \nabla \mathbf{v} : (\nabla \Delta\mathbf{u} \cdot \mathbf{S}(\mathbf{u})) \, d\Omega $$ {#eq-full-jacobian} Итерационный процесс продолжается до тех пор, пока невязка не станет меньше заданного допуска. ## Практическая задача ### Постановка задачи и импорт геометрии Мы рассмотрим стальной кронштейн (шарнир), левый торец которого жестко заделан, а на правый наложено кинематическое граничное условие — поворот на заданный угол и смещение. Вместо ручного построения геометрии мы импортируем готовую модель из файла САПР `model.step`. ```{julia} #| echo: true #| output: false #| lst-cap: "Импорт CAD геометрии и построение сетки" #| lst-label: lst-gmsh-cad gmsh.initialize() gmsh.option.setNumber("General.Terminal", 0) gmsh.model.add("model") # Импорт геометрии из STEP файла path = "model.step" |> tcf gmsh.model.occ.importShapes(path) gmsh.model.occ.synchronize() # Получение 3D объектов и их габаритов для автоматической разметки d3, = gmsh.model.getEntities(3) xmin, ymin, zmin, xmax, ymax, zmax = gmsh.model.occ.getBoundingBox(d3...) EPS = 1e-5 # Поиск поверхностей, лежащих на левом и правом краях по координате Z left_bnd = gmsh.model.getEntitiesInBoundingBox(-Inf, -Inf, zmin - EPS, +Inf, +Inf, zmin + EPS) right_bnd = gmsh.model.getEntitiesInBoundingBox(-Inf, -Inf, zmax - EPS, +Inf, +Inf, zmax + EPS) # Назначение физических групп (0, 1, 2 - точки, линии, поверхности) gmsh.model.addPhysicalGroup(3, [d3[2]], -1, "domain") for dim in (0,1,2) gmsh.model.addPhysicalGroup(dim, [tag[2] for tag in left_bnd if tag[1] == dim], -1, "fixed") gmsh.model.addPhysicalGroup(dim, [tag[2] for tag in right_bnd if tag[1] == dim], -1, "loaded") end gmsh.model.occ.synchronize() # Настройки генератора сетки gmsh.option.setNumber("Mesh.Algorithm", 4) gmsh.option.setNumber("Mesh.MeshSizeMax", 5.0) gmsh.model.mesh.generate(3) msh_file = "hinge_model.msh" |> tcf gmsh.write(msh_file) gmsh.finalize() ``` ### Дискретизация и математическая модель в Gridap Зададим свойства материала (сталь, $E = 2.1\cdot 10^5$ МПа, $\nu = 0.3$) и переведем формулы тензорной кинематики в синтаксис Gridap. Мы будем активно использовать ленивые операторы композиции `∘` и внутреннего произведения `⊙`. ```{julia} #| echo: true #| output: false #| lst-cap: "Определение нелинейной кинематики и функционалов" #| lst-label: lst-kinematics model = GmshDiscreteModel(msh_file) E_mod = 2.1e5 # Масштабированный модуль ν = 0.3 λ = (E_mod * ν) / ((1 + ν) * (1 - 2 * ν)) μ = E_mod / (2 * (1 + ν)) # Симметричная часть тензора sym(v) = 0.5 * (v + v') # Градиент деформации F F(∇u) = one(∇u) + ∇u' # Тензор деформаций Грина-Лагранжа E Ε(∇u) = 0.5 * (F(∇u)' ⋅ F(∇u) - one(F(∇u))) # Вариация тензора деформаций dE dΕ(∇du, ∇u) = sym(∇du ⋅ F(∇u)) # Второй тензор напряжений Пиолы-Кирхгофа S (модель Сен-Венана - Кирхгофа) S(∇u) = λ * tr(Ε(∇u)) * one(Ε(∇u)) + 2 * μ * Ε(∇u) # Вариация тензора напряжений dS dS(∇du, ∇u) = λ * tr(dΕ(∇du, ∇u)) * one(dΕ(∇du, ∇u)) + 2 * μ * dΕ(∇du, ∇u) degree = 2 Ω = Triangulation(model) dΩ = Measure(Ω, degree) # Невязка (residual) res(u, v) = ∫( (dΕ ∘ (∇(v), ∇(u))) ⊙ (S ∘ ∇(u)) ) * dΩ # Якобиан: материальная и геометрическая части jac_mat(u, du, v) = ∫( (dΕ ∘ (∇(v), ∇(u))) ⊙ (dS ∘ (∇(du), ∇(u))) ) * dΩ jac_geo(u, du, v) = ∫( ∇(v) ⊙ ( (S ∘ ∇(u)) ⋅ ∇(du) ) ) * dΩ jac(u, du, v) = jac_mat(u, du, v) + jac_geo(u, du, v) ``` ### Настройка решателя и инкрементальное нагружение В нелинейных задачах полное приложение нагрузки за один шаг часто приводит к расходимости метода Ньютона-Рафсона. Чтобы обойти эту проблему, мы разобьем кинематическое нагружение на несколько инкрементальных шагов. ```{julia} #| echo: true #| output: false #| lst-cap: "Инкрементальное решение задачи с использованием метода Ньютона" #| lst-label: lst-solver using LineSearches: BackTracking reffe = ReferenceFE(lagrangian, VectorValue{3,Float64}, 1) V = TestFESpace(model, reffe, conformity=:H1, dirichlet_tags=["fixed", "loaded"]) # Настройка нелинейного решателя nls = NLSolver(show_trace=false, method=:newton, linesearch=BackTracking(), iterations=20) solver = FESolver(nls) function solve_step(x0, angle, cache) # Нулевое перемещение на "fixed" g0 = VectorValue(0.0, 0.0, 0.0) # Заданное кинематическое вращение и смещение на "loaded" g1(x) = VectorValue( x[1] * cos(angle) + x[3] * sin(angle) - x[1], 0.0, -x[1] * sin(angle) + x[3] * cos(angle) - x[3] ) U = TrialFESpace(V, [g0, g1]) op = FEOperator(res, jac, U, V) uh_step = FEFunction(U, x0) uh_step, cache = solve!(uh_step, solver, op, cache) return uh_step, get_free_dof_values(uh_step), cache end disp_max = deg2rad(45) # Поворот на 45 градусов nsteps = 5 x_vals = zeros(Float64, num_free_dofs(V)) cache = nothing uh_final = nothing for step in 1:nsteps current_angle = step * disp_max / nsteps uh_final, x_vals, cache = solve_step(x_vals, current_angle, cache) end ``` ### Визуализация результатов с помощью Makie Визуализируем норму вектора перемещений на расчетной модели. ```{julia} #| echo: true #| output: true #| label: fig-displacement #| fig-cap: "Распределение нормы перемещений в кронштейне при больших деформациях" #| lst-cap: "Визуализация поля перемещений" #| lst-label: lst-vis-disps-2d fig = Figure(size = (600, 600)) ax = Axis3(fig[1,1], aspect = :data) # Вычисляем норму перемещений u_mag = sqrt ∘ (uh_final ⋅ uh_final) u_max = maximum(u_mag.(Ω.grid.node_coordinates)); ∂Ω = BoundaryTriangulation(model) # Отрисовка геометрии с наложением цветовой карты plt = plot!(ax, ∂Ω, u_mag, colormap=cgrad(:turbo, 10, categorical=true), colorrange=(0,0.99*u_max)) wireframe!(ax, ∂Ω, color=:black, linewidth=1.2) Colorbar(fig[1, 2], plt, label = "мм") writevtk(Ω, tcf("disps"); cellfields=Dict("uh"=>uh_final, "u_mag"=>u_mag)) # save("fig.png", fig) fig ``` ### Анализ и интерпретация результатов Использование геометрически нелинейной постановки позволило корректно рассчитать поведение кронштейна при больших поворотах (до 45 градусов). Ключевые выводы: Линейная модель (малые деформации $\boldsymbol{\varepsilon} = \frac{1}{2}(\nabla \mathbf{u} + \nabla \mathbf{u}^T)$) в такой задаче привела бы к возникновению огромных фиктивных напряжений и нефизичному увеличению объема детали в зоне вращения. Использование тензора Грина-Лагранжа $\mathbf{E}$ отфильтровывает повороты как жесткого целого, позволяя материалу деформироваться корректно. Инкрементальный подход (разбиение нагрузки на 5 шагов) — критически важная техника. Матрица Якоби зависит от текущих перемещений $\mathbf{u}$, и при попытке приложить полный угол поворота за одну итерацию алгоритм Ньютона-Рафсона ушел бы в бесконечность из-за плохой начальной догадки. ## Контрольные вопросы и задания ### Теоретические вопросы В чем физический смысл разницы между градиентом перемещений $\nabla \mathbf{u}$ и тензором Грина-Лагранжа $\mathbf{E}$? Почему модель материала Сен-Венана — Кирхгофа нельзя применять для эластомеров (резин) при сильном растяжении, несмотря на то, что она допускает большие перемещения? Каков физический и математический смысл геометрической матрицы жесткости в уравнении Якобиана? ### Практические задания 1. Измените количество шагов инкрементального нагружения nsteps до 1. Изучите вывод решателя. Произойдет ли сходимость? 2. Измените угол поворота disp_max на 60 и 90 градусов. Как изменяется требуемое количество итераций на шаг? 3. (\*) Реализуйте и визуализируйте расчет эквивалентных напряжений по Мизесу (через тензор Коши, который получается путем пуш-форварда тензора Пиолы-Кирхгофа: $\boldsymbol{\sigma} = J^{-1} \mathbf{F} \mathbf{S} \mathbf{F}^T$).