using Arpack, SparseArrays, LinearAlgebra, CCMechBook
# Геометрические параметры (генерация скрыта в приложении)
w, h, l = 0.1, 0.03, 0.7
model = generate_link_model(w, h, l) # 3D модель тяги
# Свойства материала (Алюминий)
E, ν, ρ = 70e9, 0.33, 2700.0
λ = (E*ν)/((1+ν)*(1-2*ν)); μ = E/(2*(1+ν))
σ(ε) = λ*tr(ε)*one(ε) + 2*μ*ε
# Пространства БЕЗ ГУ (полный контроль над DOF)
reffe = ReferenceFE(lagrangian, VectorValue{3,Float64}, 1)
V = TestFESpace(model, reffe, conformity=:H1)
U = TrialFESpace(V)
Ω = Triangulation(model); dΩ = Measure(Ω, 2)
blf_stiff(u, v) = ∫(ε(v) ⊙ (σ ∘ ε(u))) * dΩ
blf_mass(u, v) = ρ * ∫(v ⋅ u)dΩ
mK = assemble_matrix(blf_stiff, V, U)
mM = assemble_matrix(blf_mass, V, U)4 Модальный анализ и кинематические связи (RBE2)
4.1 Аннотация и цели обучения
Физическая задача: Расчет собственных частот и форм колебаний 3D-модели тяги (звена механизма) с двумя крепежными отверстиями.
Освоенные численные методы: Сборка матриц без наложения граничных условий, имплементация метода конденсации для создания абсолютно жестких связей (RBE2), применение условий Дирихле через редуцирование системы (исключение блоков), визуализация деформированного состояния.
Предварительные знания: 3D теория упругости, понятие билинейных форм, основы динамики механических систем
4.2 Физическая теория: Динамика сплошной среды
4.2.1 Сильная постановка
В отличие от задач статики, где рассматривается равновесие внутренних и внешних сил, в задачах динамики мы должны учитывать силы инерции и силы сопротивления (демпфирования). Уравнение движения сплошной среды (эластодинамики) имеет вид:
\[ \nabla \cdot \boldsymbol{\sigma} + \mathbf{f} = \rho \mathbf{\ddot{u}} + \xi \mathbf{\dot{u}} \]
где \(\boldsymbol{\sigma}\) — тензор напряжений, \(\mathbf{f}\) — объемные силы, \(\rho\) — плотность, \(\xi\) — коэффициент вязкого трения, \(\mathbf{\ddot{u}}\) и \(\mathbf{\dot{u}}\) — ускорение и скорость соответственно.
4.2.2 Слабая постановка
Умножая уравнение на произвольную пробную функцию (виртуальное перемещение) \(\mathbf{v}\) и интегрируя по объему тела \(\Omega\), с применением формулы интегрирования по частям мы получаем слабую форму:
\[ \int_{\Omega} \rho (\mathbf{v} \cdot \mathbf{\ddot{u}}) \, d\Omega + \int_{\Omega} \xi (\mathbf{v} \cdot \mathbf{\dot{u}}) \, d\Omega + \int_{\Omega} \boldsymbol{\varepsilon}(\mathbf{v}) : \boldsymbol{\sigma}(\mathbf{u}) \, d\Omega = \int_{\Omega} \mathbf{v} \cdot \mathbf{f} \, d\Omega \tag{4.1}\]
4.2.3 Матричная запись
После дискретизации области конечными элементами эта интегральная форма превращается в систему обыкновенных дифференциальных уравнений второго порядка:
\[ \mathbf{M}\mathbf{\ddot{u}} + \mathbf{C}\mathbf{\dot{u}} + \mathbf{K}\mathbf{u} = \mathbf{F} \tag{4.2}\]
где \(\mathbf{K}\) — матрица жесткости, \(\mathbf{M}\) — матрица масс и \(\mathbf{C}\) — матрица демпфирования. Матрица жесткости формируется из билинейной формы \(a(\mathbf{u}, \mathbf{v}) = \int_{\Omega} \boldsymbol{\varepsilon}(\mathbf{v}) : \boldsymbol{\sigma}(\mathbf{u}) \, d\Omega\) а матрица масс - из \(m(\mathbf{u}, \mathbf{v}) = \int_{\Omega} \rho (\mathbf{v} \cdot \mathbf{u}) \, d\Omega\). С матрицей демпфирования дела обстоят несколько сложнее: хотя в уравнении (Уравнение 4.1) присутствует билинейная форма для демпфирования \(с(\mathbf{u}, \mathbf{v}) = \int_{\Omega} \xi (\mathbf{v} \cdot \mathbf{u}) \, d\Omega\), однако на практике матрица \(\mathbf{C}\) редко собирается напрямую из нее, так как микроскопические механизмы трения сложны для точного описания. Вместо этого используется модель рэлеевского демпфирования, где матрица \(\mathbf{C}\) генерируется как линейная комбинация: \(\mathbf{C} = \alpha\mathbf{M} + \beta\mathbf{K}\), а коэффициенты \(\alpha\) и \(\beta\) определяются экспериментально или из практических соображений: коэффициент \(\alpha\) отвечает за демпфирование, возникающее в следствие взаимодействия с внешней средой (воздуха или жидкости) и влияющее больше на низкие частоты, а коэффициент \(\beta\) отвечает за демпфирование связанное с внутренним трением материала и влияющее больше на высокие частоты колебаний.
4.3 Математическая теория: Проблема собственных значений
При анализе собственных частот мы исследуем свободные колебания системы без учета внешних сил (\(\mathbf{F} = \mathbf{0}\)) и демпфирования (\(\mathbf{C} = \mathbf{0}\)). Тогда уравнение (Уравнение 4.2) примет вид:
\[ \mathbf{M}\mathbf{\ddot{u}} + \mathbf{K}\mathbf{u} = 0 \tag{4.3}\]
Известно, что такое уравнение имеет аналитическое решение в виде гармонических колебаний: \(\mathbf{u}(t) = \boldsymbol{\phi} \sin(\omega t)\). Подставляя это решение (и его вторую производную \(\mathbf{\ddot{u}} = -\omega^2 \boldsymbol{\phi} \sin(\omega t)\)) в уравнение движения, получаем:
\[ -\omega^2 \mathbf{M} \boldsymbol{\phi} + \mathbf{K} \boldsymbol{\phi} = \mathbf{0} \quad \implies \quad \mathbf{K} \boldsymbol{\phi} = \lambda \mathbf{M} \boldsymbol{\phi} \tag{4.4}\]
где \(\lambda = \omega^2\) — собственное значение, \(\boldsymbol{\phi}\) — собственный вектор (форма колебаний), а \(\omega = 2\pi f\) — круговая частота.
Равенство (Уравнение 4.4) носит название обобщенной проблемы собственных значений и в общем случае имеет \(N\) независимых решений \((\lambda_i,{ }\phi_i)\). Данные решения обладают следующими свойствами: - Собственные формы \(\boldsymbol{\phi}_i\) и \(\boldsymbol{\phi}_j\) ортогональны относительно матрицы масс: \(\boldsymbol{\phi}_i^T \mathbf{M} \boldsymbol{\phi}_j = 0\) при \(i \neq j\). Обычно их нормируют так, чтобы \[ \boldsymbol{\Phi}^T \mathbf{M} \boldsymbol{\Phi} =\mathbf{I}. \tag{4.5}\]
Однако может встречатся обычная нормировка, когда сумму квадратов всех компонентов равна единице, т.е. \(\boldsymbol{\Phi}^T \boldsymbol{\Phi} = \mathbf{I}\). Такой вид нормирования зачастую используется для визуализации собственных форм колебаний.
Кроме того, формы ортогональны относительно матрицы жесткости: \[ \boldsymbol{\Phi}^T \mathbf{K} \boldsymbol{\Phi} = \begin{bmatrix} \omega_1^2&0&\cdots&0\\ 0&\omega_2^2&\cdots&0\\ \vdots&\vdots&\ddots&\vdots\\ 0&0&\cdots&\omega_n^2\\ \end{bmatrix}. \]
Данный факт крайне удобен для анализа вынужденных колебаний, так как позволяет разложить решение на независимые моды и заменить систему уравнений на \(N\) независимых уравнений второго порядка, коэффициент жесткости в которых равен \(\omega_i^2\) а масса равна единице:
\[ \ddot{q}_i + \omega_i^2 q_i = F_i(t), \quad i = 1,2,\ldots,N{,} \]
где \(q_i\) — обобщенная координата, а \(F_i(t)\) – модальная сила, вычисляемая как проекция внешней силы на собственную форму \(\boldsymbol{\phi}_i\).
4.3.1 Численное решение
Для реальных КЭ-моделей размерность матриц составляет \(10^5-10^7\), поэтому искать все \(N\) собственных значений невозможно. Нас интересуют только первые \(n\) низших частот, т.к. низкие частоты имеют наибольший вклад в движение системы. Для поиска низших сосбтвенных значений и соответствующих им векторов используются итерационные алгоритмы (алгоритм Арнольди или Ланцоша). Наиболее развитые реализации эти алгоритмов можно найти в библиотеке Arpack.
Однако даже с применением итерационных методов поиск наименьших по амплитуде собственных значений для пучка \((\mathbf{K}{,}\mathbf{M})\) может быть проблематичным и алгоритм может не сойтись. Поэтому на практике часто ищут наибольшие собственные значения обратного пучка \((\mathbf{M}{,}\mathbf{K})\):
\[ \mathbf{M} \boldsymbol{\phi} = \lambda \mathbf{K}\boldsymbol{\phi}, \quad \text{где } \lambda = \frac{1}{\omega^2} \]
Найденные наибольшие \(\lambda\) соответствует искомым низшим частотам \(f = \frac{1}{2\pi \sqrt{\lambda}}\). Собственные вектора ф этом случае будут отнормированы относительно матрицы \(\mathbf{K}\), и при необходимости их можно реортонормировать относительно \(\mathbf{M}\). Поскольку вектора уже ортогональны, то для этого достаточно их отмасштабировать. Для этого необходимо вычислить эффективную модальную массу для каждого \(i\)-го вектора по формуле \(m_i = \phi_i^T \mathbf{M} \phi_i\). Затем каждый исходный вектор нормируется путем деления на корень из полученного значения:\[\psi_i = \frac{\phi_i}{\sqrt{m_i}}\]
Отдельный интерес представляет задача поиска собственных значений в случае вырожденных систем, когда матрица \(\mathbf{K}\) имеет нулевую собственную частоту. В этом случае алгоритм Арнольди может не сойтись, и необходимо использовать метод обратных итераций (в англоязычной литературе данный метод носит название shift-invert). Суть данного метода заключается в следующем: вместо поиска собственных значений пучка \((\mathbf{M}{,}\mathbf{K})\) ищутся собственные значения сдвинутого пучка \((\mathbf{M}{,} \mathbf{B})\), где \(\mathbf{B}=\mathbf{K}-\sigma \mathbf{M}\), a \(\sigma\) – положительное число, близкое к нулю. Собственные значения \(\mu\) пучка \((\mathbf{M}{,} \mathbf{B})\) связаны c собственными значениями \(\lambda\) исходного пучка соотношением: \[ \mu = \frac{1}{\lambda - \sigma}. \]
После нахождения собственных значений сдвинутого пучка, они преобразуются обратно в собственные значения исходного пучка с помощью формулы \(\lambda = \frac{1}{\mu} + \sigma\).
В библиотеке Arpack.jl данный метод реализован через ключевой аргумент sigma функции eigs. При этом, функция возращает уже восстановленные собственные значения \(\lambda\) и собственные вектора \(\boldsymbol{\phi}\) исходного пучка \((\mathbf{K}{,}\mathbf{M})\). Если параметр sigma не указан, то используется стандартный алгоритм Арнольди.
4.4 Кинематические связи
В инженерной практике интерфейсные области, используемые для соединения деталей (отверстия под болты, подшипники) часто моделируют с помощью кинематических элементов. В литературе их можно встретить под названием Rigid Body Elements (RBE), введенные в инженерном ПО NASTRAN, которые с тех пор используются как полноценные инженерные термины. При этом различают абсолютно жесткий элемент (RBE2) и элемент распределения нагрузок (RBE3), в общем случае, не добавляющий дополнительной жесткости в систему. Для каждого из этих элементов характерно наличие управляющих и зависимых узлов/степеней свободы. Для обозначения управляющего узла в литературе часто используют термин master node, а для зависимого узла — slave node. В случае RBE2 зависимые узлы жестко связаны с управляющим узлом, что позволяет передавать как поступательные перемещения, так и вращения.
Математически это означает, что перемещение любого зависимого узла \(\mathbf{u}_s\) определяется перемещением и поворотом управляющего узла \(\mathbf{u}_m\):
\[ \mathbf{u}_s = \mathbf{u}_m + \boldsymbol{\theta}_m \times \mathbf{r} \]
где \(\mathbf{r}\) — вектор от управляющего узла к зависимому узлу. В общем случае, когда у узлов также присутствуют поворотные степени свободы (например, узлы балочных или оболочечных элементов), так же добавляются уравнения связывающие углы управляющего и зависимых узлов: \[ \mathbf{\theta}_s = \mathbf{\theta}_m \]
Для реализации этого в МКЭ строится матрица трансформации (проекции) \(\mathbf{C}\), выражающая полный вектор степеней свободы \(\mathbf{u}\) через сокращенный вектор независимых степеней свободы \(\tilde{\mathbf{u}}\):
\[ \mathbf{u} = \mathbf{C} \tilde{\mathbf{u}} \]
Подставляя это в уравнение колебаний Уравнение 4.3 мы получим уравнение уравнение вида: \[ \mathbf{MC}\ddot{\tilde{\mathbf{u}}} + \mathbf{KC}\tilde{\mathbf{u}} = 0. \tag{4.6}\]
После этого мы можем применить проекцию Галеркина и потребовать ортогональности ошибки к матрице \(mathbf{C}\), для чего умножим уравнение Уравнение 4.6 слева на \(\mathbf{C}^T\). В результате получим редуцированную систему:
\[ (\mathbf{C}^T \mathbf{K} \mathbf{C}) \tilde{\boldsymbol{\phi}} = \lambda (\mathbf{C}^T \mathbf{M} \mathbf{C}) \tilde{\boldsymbol{\phi}} \]
Подробнее о методе Галеркина можно прочитать в [1].
4.5 Практический кейс: шатун с RBE2 интерфейсами
В данном примере мы реализуем “ручное” управление сборкой. Мы не задаем ГУ в пространствах Gridap, а модифицируем итоговые матрицы. Смоделируем трехмерное звено механизма из алюминиевого сплава (\(\rho = 2700\) кг/м³, \(E = 70\) ГПа, \(\nu = 0.33\)). Звено имеет длину 700 мм, высоту 80 мм и ширину 100 мм. По краям расположены отверстия под подшипники. Внешний вид шатуна представлен на Рисунок 5.1.
4.5.1 Построение пространств и базовая сборка матриц
Сформируем векторные пространства испытательных и пробных функций без использования параметра dirichlet_tags, чтобы все узлы модели имели полный набор степеней свободы – по 3 поступательных перемещения. Мы сознательно не используем стандартные граничные условия Дирихле на этапе генерации конечно-элементных пространств, поскольку необходим полный контроль над всеми степенями свободы системы. Поэтому мы сначала соберем исходные матрицы жесткости и масс, а после этого модифицируем их через матрицу связей \(\mathbf{C}\) для внедрения RBE2 элементов и наложим граничные условия через уже на эти элементы.
4.5.2 Имплементация связей RBE2 (Конденсация матриц)
Мы связываем узлы, принадлежащие цилиндрической поверхности отверстия с мастер-узлами, лежащими на оси отверстия. Это уменьшает количество степеней свободы и корректно передает моменты. В русской литературе этот подход часто называют методом конденсации (исключение зависимых узлов).
left_nodes = CCMechBook.get_nodes_by_tag(model, "left_hole")
right_nodes = CCMechBook.get_nodes_by_tag(model, "right_hole")
master_coords = [VectorValue(0.0, 0.0, 0.0), VectorValue(l, 0.0, 0.0)]
slave_nodes = [left_nodes, right_nodes]
n_glob = num_free_dofs(V)
dofs_slave_all = Int64[]
for node in [right_nodes; left_nodes]
append!(dofs_slave_all, V.metadata.node_and_comp_to_dof[node])
end
dofs_indep = setdiff(1:n_glob, dofs_slave_all)
# Матрица C: [N_glob x N_reduced]
C_rbe = spzeros(n_glob, length(dofs_indep) + 12)
C_rbe[dofs_indep, 1:length(dofs_indep)] = I(length(dofs_indep))
cur_col = length(dofs_indep) + 1
node_coords = get_node_coordinates(Ω)
for i in 1:2
for node in slave_nodes[i]
r = node_coords[node] - master_coords[i]
s_dofs = V.metadata.node_and_comp_to_dof[node]
C_rbe[Vector(s_dofs), cur_col:cur_col+2] = I(3)
C_rbe[Vector(s_dofs), cur_col+3:cur_col+5] = -CCMechBook.skew_symmetric(r)
end
global cur_col += 6
end
mK_cond = C_rbe' * mK * C_rbe
mM_cond = C_rbe' * mM * C_rbe4.5.3 Наложение граничных условий и решение
Для консольной заделки мы исключим из матриц строки и столбцы, соответствующие степеням свободы левого мастер-узла.
n_cond = size(C_rbe, 2)
left_rbe_dofs = (n_cond - 11):(n_cond - 6) # DOF левого отверстия
# Исключаем блоки, соответствующие заделке (все 6 DOF)
free_dofs = setdiff(1:n_cond, left_rbe_dofs)
mK_final = mK_cond[free_dofs, free_dofs]
mM_final = mM_cond[free_dofs, free_dofs]Для решения задачи \(\mathbf{K}\boldsymbol{\phi} = \omega^2 \mathbf{M}\boldsymbol{\phi}\) используем пакет Arpack.jl, который реализует итерационный метод неявно перезапущенного алгоритма Арнольди.
num_modes = 10
λ, Φ_red = eigs(mM_final, mK_final; nev=num_modes, which=:LM)
# Восстановление векторов до полной размерности
Φ_full = zeros(n_cond, num_modes)
Φ_full[free_dofs, :] = Φ_red
Φ_final = C_rbe * Φ_full
# Экспорт в VTK
cellfields = Dict{String,Any}()
for i in 1:num_modes
f = round(sqrt(1 / abs(λ[i])) / (2π), digits=2)
cellfields["mode_$(i)_f=$(f)Hz"] = FEFunction(U, real.(Φ_final[:, i]))
end
writevtk(Ω, tcf("eigen_modes"); cellfields)4.5.4 Визуализация и анализ результатов
На Рисунок 4.2 представлены первые 4 собственные формы колебаний: первая изгибная форма колебаний из рабочей плоскости шатуна, первая изгибная форма в рабочей плоскости, вторая изгибная форма из рабочей плоскости и первая крутильная форма колебаний.
По представленным формам колебаний хорошо видно условия закрепления: на всех формах можно наблюдать, что перемещения левого конца балки нулевые; по первой и третьей формам видно, что отсутствует поворот закрепленного конца балка в плоскости \(XZ\), равно как по второй форме видно отсутствие поворота закрепленного конца балки в плоскости \(XY\); на четвертой форме отчетливо видно нулевой угол поворота закрепленного конца в плоскости \(YZ\). Кроме того, на третьей собственной форме можно наблюдать узел колебаний — область балки, в которой отсутствуют перемещения, хотя там нет никакого закрепления. Эту область отчетливо видно по концентрическим окружностям на цветовой карте перемещений. Так же обратите внимание на центральную ось балки при кручении — она остается неподвижной на всей длине, что хорошо видно по темносинему цвету на карте перемещений.
4.6 Контрольные вопросы и задания
4.6.1 Теоретические вопросы
Почему для свободного тела в плоской (2D) постановке первые 3 собственные частоты равны нулю, а в пространственной (3D) — первые 6?
Почему в МКЭ для динамики выгодно искать собственные значения пучка \(\mathbf{M}{,}\mathbf{K}\) соответствующего матрице \(\mathbf{K}^{-1}\mathbf{M}\), а не наоборот?
Как изменится собственная частота первой изгибной формы, если материал звена заменить с алюминия на сталь (модуль упругости выше в 3 раза, плотность выше в 2.9 раза)? Подтвердите ответ теоретической формулой \(\omega \propto \sqrt{E/\rho}\).
4.6.2 Практические задания
Сравните значения полученных частот со значениями полученными аналитически для балки прямоугольного поперечного сечения и вычислите погрешность. Объясните чем обсуловлена данная погрешность и какая модель более точная?
Освободите вращательную степень свободы (вокруг оси \(Z\)) левого мастер-узла и закрепите аналогичным образом правый конец балки. Как изменятся первые 4 формы колебаний?
Используя исходные условия закрепления, Измените плотность материала и проверьте, соблюдается ли зависимость \(f \propto 1/\sqrt{\rho}\).
(*) Используя исходные условия закрепления, к правому отверстию присоедените сосредоточеную массу. Постройте график зависимости частоты первой изгибной моды от массы груза \(M_{load}\) (варьируйте от 1 до 20 кг). Сравните с аналитическим расчетом для балки с грузом на конце.
4.7 Приложение
function generate_link_model(w, h, l)
EL_ORDER = 1
hole_diam = 50mm
gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 1)
gmsh.option.setNumber("Mesh.ElementOrder", EL_ORDER)
gmsh.model.add("link_mesh")
lc = 0.025
p1 = gmsh.model.occ.add_point(0.0, -w / 2, -h / 2, lc)
p2 = gmsh.model.occ.add_point(l, -w / 2, -h / 2, lc)
p3 = gmsh.model.occ.add_point(l, w / 2, -h / 2, lc)
p4 = gmsh.model.occ.add_point(0.0, w / 2, -h / 2, lc)
p5 = gmsh.model.occ.add_point(0.0, 0.0, -h / 2, lc)
p6 = gmsh.model.occ.add_point(l, 0.0, -h / 2, lc)
p7 = gmsh.model.occ.add_point(hole_diam / 2, 0.0, -h / 2, lc)
p8 = gmsh.model.occ.add_point(-hole_diam / 2, 0.0, -h / 2, lc)
p9 = gmsh.model.occ.add_point(l - hole_diam / 2, 0.0, -h / 2, lc)
p10 = gmsh.model.occ.add_point(l + hole_diam / 2, 0.0, -h / 2, lc)
c1 = gmsh.model.occ.add_circle_arc(p1, p5, p4, -1, true)
c2 = gmsh.model.occ.add_circle_arc(p3, p6, p2, -1, true)
c3 = gmsh.model.occ.add_circle_arc(p8, p5, p7, -1, true)
c4 = gmsh.model.occ.add_circle_arc(p7, p5, p8, -1, true)
c5 = gmsh.model.occ.add_circle_arc(p9, p6, p10, -1, true)
c6 = gmsh.model.occ.add_circle_arc(p10, p6, p9, -1, true)
l1 = gmsh.model.occ.add_line(p1, p2)
l2 = gmsh.model.occ.add_line(p3, p4)
curve_loop = gmsh.model.occ.add_curve_loop([-c1, l1, -c2, l2])
hole1 = gmsh.model.occ.add_curve_loop([c3, c4])
hole2 = gmsh.model.occ.add_curve_loop([c5, c6])
surface = gmsh.model.occ.add_plane_surface([curve_loop, hole1, hole2])
volume = gmsh.model.occ.extrude([(2, surface)], 0, 0, h)
gmsh.model.occ.synchronize()
EPS = 1e-6
hole1_bb = (
-hole_diam / 2 - EPS,
-hole_diam / 2 - EPS,
-h / 2 - EPS,
hole_diam / 2 + EPS,
hole_diam / 2 + EPS,
h / 2 + EPS)
hole2_bb = (l - hole_diam / 2 - EPS,
-hole_diam / 2 - EPS,
-h / 2 - EPS,
l + hole_diam / 2 + EPS,
hole_diam / 2 + EPS,
h / 2 + EPS,
)
#
hole1_dt = gmsh.model.occ.get_entities_in_bounding_box(hole1_bb...)
hole2_dt = gmsh.model.occ.get_entities_in_bounding_box(hole2_bb...)
#
setdiff!(hole1_dt, [(0, p5)])
setdiff!(hole2_dt, [(0, p6)])
add_boundary_label(hole1_dt, "left_hole")
add_boundary_label(hole2_dt, "right_hole")
#
dim, tag = gmsh.model.occ.get_entities(3)[1]
#
gmsh.model.add_physical_group(dim, [tag], -1, "domain")
gmsh.model.occ.synchronize()
gmsh.model.mesh.generate(3)
gmsh.model.occ.synchronize()
# gmsh.fltk.run()
#
gmsh.write("link_mesh.msh" |> tcf)
gmsh.finalize()
model = GmshDiscreteModelNoVerbose("link_mesh.msh" |> tcf)
# writevtk(model, "link_model" |> tcf)
return model
end