Engee 文档
Notebook

部件稳态热场计算：Julia、Gmsh、Gridap、VTK 和 Makie

让我们应用一系列工具，在 Engee 中进行设计、设置模型参数、应用计算网格、构建要通过有限元法解决的问题，并最终计算复杂部件的传热。

简介

对于需要精确求解偏方程的问题（如传热或变形体力学），专业工具（如 Gridap 或 FEniCS）能提供更精确、更灵活的结果。它们可以处理任意网格和复杂的边界条件，这对科学计算至关重要。

In [ ]: 
Pkg.add(["Gridap", "GridapGmsh", "CairoMakie", "WGLMakie", "WriteVTK", "ReadVTK", "GeometryBasics"])

在本例中，我们使用了一个库，该库的安装会自动将 Gmsh 软件包添加到系统中。

In [ ]: 
using GridapGmsh

让我们从该库的资源库中提取一个简单的 模型示例，并稍作修改：

-boundary 1 区域将结合较宽孔的内表面（右图）

  • 区域boundary 2 - 带有突起的整个孔内表面（左图）

screenshot_from_2025_05_11_12_50_03.png

这种参数化模型的创建可以在文本编辑器中完成，例如通过在 https://gmsh.info/ 或众多 CAD 系统（包括开源 CAD 系统）中对结果进行可视化：OpenCASCADE、OpenSCAD 和 FreeCAD。

让我们以 .geo 格式（Gmsh 中几何体的基本格式）加载示例并创建网格：

In [ ]: 
GridapGmsh.gmsh.initialize()
GridapGmsh.gmsh.clear()
GridapGmsh.gmsh.parser.parse("$(@__DIR__)/demo.geo");
GridapGmsh.gmsh.model.geo.synchronize();
GridapGmsh.gmsh.option.setNumber("Mesh.CharacteristicLengthMin", 0.01)  # минимальный размер
GridapGmsh.gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 0.05)  # максимальный размер
GridapGmsh.gmsh.model.mesh.generate(3);
GridapGmsh.gmsh.write("demo.msh");

让我们输出图案的模型，这是一组四面体（Makie 不允许使用四面体，因此我们手动将其转换为三角形）。图形的颜色与点的 X 坐标相对应。

In [ ]: 
using Gridap, GridapGmsh, GeometryBasics

import CairoMakie: mesh, NoShading
# import WGLMakie: mesh, NoShading # На замену CairoMakie - для интерактивных графиков

model = GmshDiscreteModel("demo.msh");

grid = get_grid(get_triangulation(model))
points = reinterpret(reshape, Float64, grid.node_coordinates)
cells = hcat(grid.cell_node_ids...)

tetrahedra = [ (c[1], c[2], c[3], c[4]) for c in eachcol(cells) ]
all_faces = [ [TriangleFace{Int}(ta, tb, tc),TriangleFace{Int}(ta, tb, td),TriangleFace{Int}(ta, tc, td),TriangleFace{Int}(tb, tc, td)] for (ta, tb, tc, td) in tetrahedra ]
all_faces = reshape(reinterpret(Int, vcat(all_faces...)), 3, :)
colors = points[1,:]

mesh(points', all_faces', color=colors, shading=NoShading)

Info    : Reading 'demo.msh'...
Info    : 188 entities
Info    : 3302 nodes
Info    : 12065 elements
Info    : Done reading 'demo.msh'

Warning : Gmsh has aleady been initialized
Out[0]:
No description has been provided for this image

数值问题表述

我们使用基本示例中的 GridapGmsh 代码，用有限元法（FEM）求解弱（变分法）形式的带 Dirichlet 边界条件的拉普拉斯方程。

拉普拉斯方程和边界条件： 我们在区域$\Omega$ 中搜索满足方程$\Delta u = 0$ 的函数$u(x)$ ，边界条件$u = 0$ at$\Gamma_1$,$u = 1$ at$\Gamma_2$ ，其中$\Omega$ 是网格定义的区域demo.msh,$\Gamma_1$ 和$\Gamma_2$ 分别是标记为 "boundary1 "和 "boundary2 "的边界。

Slack form（变分法问题表述） 是现代有限元的一种关键方法，它放宽了对解法平滑性的要求，而且可以直接给出边界条件，而无需对其进行近似。当然，它也是问题离散化的主要方法，它允许我们在基函数中指定问题，并将其作为线性系统求解$Au = b$ 。因此，我们的问题在变分公式中被转化为$a(u, v) = \int_{\Omega} \nabla v \cdot \nabla u \, d\Omega = 0 \quad \forall v \in V_0$ ，其中$\nabla$ 是梯度，$\Omega$ 是积分区域，$V_0$ 是在边界变为零的试验空间。

请注意∫( ∇(v)⋅∇(u) )dΩ 代码中的描述与变分公式中的拉普拉斯方程是多么接近$a(u, v) = \int_{\Omega} \nabla v \cdot \nabla u \, d\Omega = 0 \quad \forall v \in V_0$ 。

In [ ]: 
using Gridap

model = GmshDiscreteModel("demo.msh");         # Загружаем геометрию сетки из файла *.msh
order = 1
reffe = ReferenceFE(lagrangian,Float64,order)  # Линейные лагранжевы элементы 1-го порядка
V = TestFESpace(model,reffe,dirichlet_tags=["boundary1","boundary2"]) # Пробное пространство
U = TrialFESpace(V,[0,1])                      # Определение пространства решений
                                               # - на boundary1 значение решения u = 0
                                               # - на boundary2 значение u = 1
Ω = Triangulation(model)                       # Триангуляция и мера интегрирования
dΩ = Measure(Ω,2*order)
a(u,v) = ∫( ∇(v)⋅∇(u) )dΩ                      # Слабая форма уравнения Лапласа
l(v) = 0                                       # Линейная форма (правая часть)
op = AffineFEOperator(a,l,U,V)                 # Собираем линейный оператор задачи
uh = solve(op)                                 # Решение системы

writevtk(Ω,"demo",cellfields=["uh"=>uh]);

Info    : Reading 'demo.msh'...
Info    : 188 entities
Info    : 3302 nodes
Info    : 12065 elements
Info    : Done reading 'demo.msh'

我们将更进一步，除了即时可视化之外，我们还将展示如何将保存在 VTK 中的计算结果可视化。但首先，我们将展示刚刚得到的解决方案：

In [ ]: 
using GeometryBasics
nodes = GridapGmsh.get_node_coordinates(Ω)
cell_node_ids = GridapGmsh.get_cell_node_ids(Ω)
uh_values = [uh(node) for node in nodes]
points = [Point3f(x...) for x in nodes]
faces = [TriangleFace(cell_node_ids[i][1], cell_node_ids[i][2], cell_node_ids[i][3]) for i in 1:length(cell_node_ids)]
mesh(points, faces, color=uh_values, colormap=:inferno, shading=NoShading)
Out[0]:
No description has been provided for this image

可视化 VTK 文件中的解决方案

对于查看 VTK 文件，最方便的解决方案是开源程序 ParaView。但对我们来说，在 Engee 中获得可完全发现的交互式报告也很有价值。因此，我们将在 Julia 中打开该文件，并在本脚本中直接输出解决方案。

如果我们将这一部分与解决方案异步运行（可能需要很长时间），则需要连接这些库：

In [ ]: 
using WriteVTK, ReadVTK

它们将帮助我们读取 VTK 文件并将其显示在 3D 图形上。有些功能还需要额外的库，但我们现在只能使用低级命令：

  • 创建顶点坐标矩阵points
  • 获取矢量uh_values ，我们的解决方案的场uh 值将位于该处。
  • 将四面体cells 分割成一组三角形all_faces
In [ ]: 
vtk = VTKFile("demo.vtu")
points = get_points(vtk)
uh_values = get_data(get_point_data(vtk)["uh"])

cells = get_cells(vtk)

# Разобьем массив cells.connectivity на тетраедры, а затем на треугольники
using GeometryBasics
all_faces = TriangleFace{Int}[]
offset = 1
while offset <= length(cells.connectivity)
    ta, tb, tc, td = cells.connectivity[offset:offset+3]
    push!(all_faces, TriangleFace(ta, tb, tc))  # грань 1
    push!(all_faces, TriangleFace(ta, tb, td))  # грань 2
    push!(all_faces, TriangleFace(ta, tc, td))  # грань 3
    push!(all_faces, TriangleFace(tb, tc, td))  # грань 4
    offset += 4
end

# CairoMakie не принимает последний треугольник в модели, временно используем [1:end-1]
mesh(points, all_faces[1:end-1], color=uh_values, colormap=:coolwarm, shading=NoShading)
Out[0]:
No description has been provided for this image

结论

我们在四面体网格上实现了拉普拉斯方程的求解，将结果保存为 VTK 格式并将其可视化--这是物理场（热、电、机械）建模的基础。这种方法可确保实际工程问题的可重复性、准确性和可扩展性，在实际工程问题中，需要的是详细的计算而不是简化的系统近似。