详细计算静止热场：Julia，Gmsh，Gridap，VTK和Makie

让我们应用一系列工具，使Engee能够设计、参数化模型、对其应用计算网格、构建使用有限元法求解的问题，从而计算复杂零件的传热。

导言

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

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

对于这个例子，我们使用一个库，安装它会自动将Gmsh包添加到系统中。

using GridapGmsh

让我们来看一个简单的[模型示例]（https://github.com/gridap/GridapGmsh.jl/tree/master/demo ）从这个库的存储库并稍微修改它:

-地区 boundary 1 将统一较宽孔的内表面（右侧）
-地区 boundary 2 -带突起的孔的整个内表面（左侧）

您可以在文本编辑器中创建这样的参数化模型，例如，通过在https://gmsh.info /或者在许多用于固态建模的CAD系统之一中，包括开源系统：[OpenCascade]（https://dev.opencascade.org /），[OpenSCAD]（https://openscad.org /)和[FreeCAD](https://www.freecad.org /)。

让我们上传一个格式的示例。geo（Gmsh中几何的基本格式）并创建网格:

GridapGmsh.gmsh.initialize()
GridapGmsh.gmsh.clear()
GridapGmsh.gmsh.open("$(@__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轴上点的坐标。

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

设置数值问题

我们使用[基本示例]中的代码（https://github.com/gridap/GridapGmsh.jl/tree/master ）对于GridapGmsh，其中拉普拉斯方程用有限元法（FEM）在弱（变分）公式中用Dirichlet边界条件求解。

**拉普拉斯方程和边界条件：**正在搜索函数满足方程在具有边界条件上 , 上，在哪里 -网格定义的区域 demo.msh, 和 -分别带有标签"boundary1"和"boundary2"的边框。

**弱形式（问题的变分公式）**是现代FEM的关键方法，它削弱了对解平滑度的要求，可以直接设置边界条件而无需近似。而且，当然，这是问题离散化的主要方法，它允许您在基本功能中设置它并将其作为线性系统求解。 . 所以，我们在变分公式中的问题转化为，在哪里 -梯度, -整合领域, -在边界处变为零的测试空间。

请注意代码中给出的描述有多接近 ∫( ∇(v)⋅∇(u) )dΩ 变分公式中的拉普拉斯方程 .

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中的计算。但首先，让我们展示我们刚刚改进的解决方案。:

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)

从VTK文件可视化解决方案

查看VTK文件最方便的解决方案是一个开源程序[ParaView]（https://www.paraview.org /)。但对于我们来说，在Engee中完全开放的交互式报告也是有价值的。因此，我们将在Julia中打开此文件，并在此脚本中直接输出解决方案。

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

using WriteVTK, ReadVTK

他们将帮助我们读取VTK文件并将其显示在3D图形上。有一些函数有额外的库，但现在我们将使用低级命令。:

-创建矩阵 points 对于顶点的坐标
-获取向量 uh_values 字段值的位置 uh 我们的解决方案
-让我们打破四面体 cells 在一组三角形上 all_faces

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)

结论

我们在四面体网格上实现了拉普拉斯方程的解，将结果保存为VTK格式并将其可视化。

这些是物理场（热、电、机械）建模的关键阶段。这种方法确保了可重复性、准确性和随时准备扩展到需要详细计算而不是简化系统近似的实际工程任务。