在扇形区域组合 Neumann 与 Dirichlet 边界求解扩散

有曲线边界的扩散问题，难点往往不在方程本身，而在几何的分段方式和边界条件的对应关系。FiniteVolumeMethod.jl 的这个例子将一个扇形分成下直边、外圆弧和上直边：两条直边不允许通量穿过，圆弧则始终保持零值。只要这三个片段的顺序在网格与回调中一致，就能把极坐标描述的问题放到笛卡尔三角网格上求解。

本文据 SciML／FiniteVolumeMethod.jl 官方教程完整整理。原页未署个人作者；项目许可证署名为 Daniel VandenHeuvel。源文 stable 页面显示由 Documenter.jl 1.8.0 在 2025 年 1 月 1 日构建，使用 Julia 1.11.2。这是文档构建信息，不是首发日期，也没有锁定全部依赖版本。本文未运行 Julia 或复现数值实验；节点数、成功标记和结果图均明确来自官方原文。

先把方程写成网格能够使用的形式

取半径为 1、夹角为 α = π/4 的扇形，也就是 45° 区域。变量 u(r, θ, t) 满足扩散方程：

∂u/∂t = ∇²u
0 < r < 1，0 < θ < α，t > 0

θ = 0：∂u/∂θ = 0
r = 1：u = 0
θ = α：∂u/∂θ = 0

u(r, θ, 0) = 1 - r
α = π/4

两条径向直边上的条件属于 Neumann 条件。在这里，角向导数为零与法向梯度为零等价，可以写成 ∇u · n = 0，其中 n 是边界外法向。外圆弧上的 u = 0 是 Dirichlet 条件，即直接固定场值。

尽管问题用极坐标描述，代码中的微分算子采用笛卡尔坐标。利用 ∇²u = ∇·∇u，可以把方程改写为守恒形式：

∂u/∂t + ∇·q = 0
q = -∇u

负号很关键：扩散通量沿场值减小的方向流动。将其误写成正梯度，就不是同一个扩散问题了。

编辑校正：原文第三条边界条件右侧把范围写成了角度范围；按该条件的含义，θ = α 直边上的自由位置参数应是 0 < r < 1。这里按几何与实现统一表述。研究页对角度的简述也以源文明确的 π/4 为准，不解释为四分之一圆。

用三段边界建立曲边网格

三角剖分必须知道哪一段边界对应哪一种条件，不能把整个外边界塞进一个不分段的节点列表。下面为每段提供一个向量，再将它们放到总的 boundary_nodes 中：

using DelaunayTriangulation, FiniteVolumeMethod, ElasticArrays

α = π / 4
points = [(0.0, 0.0), (1.0, 0.0), (cos(α), sin(α))]
bottom_edge = [1, 2]
arc = CircularArc((1.0, 0.0), (cos(α), sin(α)), (0.0, 0.0))
upper_edge = [3, 1]
boundary_nodes = [bottom_edge, [arc], upper_edge]
tri = triangulate(points; boundary_nodes)
A = get_area(tri)
refine!(tri; max_area=1e-4A)
mesh = FVMGeometry(tri)

三个初始点依次为原点、正 x 轴上的单位半径端点，以及角度 α 对应的单位圆端点。边界沿 1 → 2 → 3 → 1 闭合：下边是 [1, 2]，圆弧从点 2 走到点 3，圆心为原点，上边是 [3, 1]。使用 CircularArc 是为了让三角剖分器认识这一段是真正的圆弧，而不是仅连接两个端点的直线。

get_area(tri) 取得区域面积。max_area=1e-4A 是 Julia 的数值字面量乘法写法，相当于 1e-4 * A，让细化后的三角形面积受到限制。对本例，区域面积是 α/2 = π/8；这是由单位半径扇形几何推得的值，用于理解尺度，并非运行报告。

原站展示的构造结果为 8179 个控制体、15990 个三角形和 24168 条边。这些数值属于原文生成环境；三角剖分算法、依赖版本及其他实现变化都可能影响具体网格，不应把这些数字作为所有环境必须完全一致的验收条件。

using CairoMakie
fig, ax, sc = triplot(tri)
fig

核对片段顺序，再绑定边界条件

网格细化会在原来两个端点之间增加许多边界节点，因此检查时不应只看最初三个顶点。使用：

get_boundary_nodes(tri)

原文的结果是含有三个元素的向量，每个元素又是一个节点索引向量。第一段从 1 到 2，第二段从 2 到 3，第三段从 3 回到 1；中间索引是细化生成的节点。这个结构确认了边界仍然保有三段及其顺序。

BoundaryConditions 接收回调元组和类型元组。两个元组中第 i 个元素，都对应网格边界的第 i 段：

lower_bc = arc_bc = upper_bc = (x, y, t, u, p) -> zero(u)
types = (Neumann, Dirichlet, Neumann)
BCs = BoundaryConditions(mesh, (lower_bc, arc_bc, upper_bc), types)

三个回调虽然都返回 zero(u)，其数学含义由对应的边界类型决定。zero(u) 还保留了与 u 相适应的数值类型。若把圆弧和直边在元组中对调，程序可能仍有语法上有效的输入，却求解了不同的物理问题。原文打印的边界对象也明确列出 (Neumann, Dirichlet, Neumann)。

指定初始场与扩散系数

初始值是 1-r。在笛卡尔坐标中，半径为 sqrt(x^2 + y^2)，因此先写成二元函数，再按网格点顺序采样。扩散系数恒为 1，积分终点取 0.1：

f = (x, y) -> 1 - sqrt(x^2 + y^2)
D = (x, y, t, u, p) -> one(u)
initial_condition = [
    f(x, y) for (x, y) in DelaunayTriangulation.each_point(tri)
]
final_time = 0.1
prob = FVMProblem(
    mesh, BCs;
    diffusion_function=D,
    initial_condition,
    final_time
)

这里沿用原文的反应扩散接口，但只提供扩散函数，没有额外反应项。one(u) 表示与当前值类型相容的 1。原站显示的问题对象包含 8179 个节点，时间区间为 (0.0, 0.1)；这里对代码做了换行排版，语义未改。

如果改用通量接口，负号仍然不能丢

原文还给出了通量形式的回调：

flux = (x, y, t, α, β, γ, p) -> (-α, -β)

这里的 α、β 和 γ 是局部线性近似 u ≈ αx + βy + γ 的系数，不是前面定义扇形角度的那个几何参数。局部梯度近似为 (α, β)，所以通量返回 (-α, -β)。

编辑说明：这段回调只是原文展示的替代建模方式。上面的 prob 已经由 diffusion_function=D 构造，仅定义名为 flux 的变量不会自动改变它。若真正切换接口，需要按所用版本的 FiniteVolumeMethod.jl 接口文档重建问题，不能声称当前示例同时运行了两种方法。

求解并比较三个时刻

原教程选择 TRBDF2 时间积分器，并用 KLUFactorization 作为线性求解方法，每隔 0.01 保存一次结果，关闭并行选项：

using OrdinaryDiffEq, LinearSolve
sol = solve(
    prob,
    TRBDF2(linsolve=KLUFactorization()),
    saveat=0.01,
    parallel=Val(false)
)

作者说，在自己的经验中，这一组合通常对这类问题表现较好。这是一条经验性说明，原文没有提供完整的跨求解器性能基准；不能据此保证它在任意规模、机器或版本上最快。

原站报告 retcode: Success，插值方式为一阶线性，并列出从 0.00 到 0.10 的 11 个保存时刻和对应的 11 组节点值。这些输出来自官方文档，不是本文执行所得。本文省略逐行展开的长数值向量，完整数值输出保留于源站。

下面按 Julia 从 1 开始的索引取第 1、6、11 帧，分别对应 t = 0、0.05 与 0.1。三个图采用相同等值线级别 0:0.01:1，便于比较：

using CairoMakie
fig = Figure(fontsize=38)
for (i, j) in zip(1:3, (1, 6, 11))
    local ax
    ax = Axis(
        fig[1, i], width=600, height=600,
        xlabel="x", ylabel="y",
        title="t = $(sol.t[j])",
        titlealign=:left
    )
    tricontourf!(
        ax, tri, sol.u[j],
        levels=0:0.01:1, colormap=:matter
    )
    tightlimits!(ax)
end
resize_to_layout!(fig)
fig

初始场在原点最高，向外圆弧降到零。两条径向边是零通量边界，外圆弧则固定为零。官方图展示了这一初始分布随扩散逐渐衰减的过程。图像可帮助检查边界是否设置反了，但不能代替误差分析、守恒检查或网格收敛研究。

复用这个例子时需要保留的检查

最值得保留的是“几何分段—边界回调—边界类型”的对应关系。改变角度、增加边界片段或更换边界条件后，应重新核对 get_boundary_nodes(tri) 与两个元组，确认每一段的物理含义。若更改网格细化阈值，应评估节点数、内存和时间成本；更密的网格并不自动证明结果收敛。

本次静态检查没有在展示代码中发现读取凭据、执行系统命令、远程下载数据、字符串求值或硬编码秘密。输入几何、初始场和参数均由代码定义。原文导入了 ElasticArrays，片段中没有直接使用它，本文为保持来源一致予以保留；依赖安装与实际兼容性尚未验证。上述检查不等于全面安全审计，也不构成“无漏洞”结论。

原文末尾的 Just the code 与前面分步代码重复。本文按解释顺序保留所有建模、检查、求解和绘图步骤，不重复粘贴第二遍；额外说明均已标明。用于科研结论之前，还需要在锁定版本的独立环境中记录收敛性与误差证据。