从局部脉冲出发求解多孔介质非线性扩散-未完纪

原文标题： Porous-Medium Equation
作者： FiniteVolumeMethod.jl 文档维护者（页面未署个人名）
来源： FiniteVolumeMethod.jl 官方教程
原文日期： 原页面未注明。页面构建信息为 Documenter.jl 1.8.0 于 2025 年 1 月 1 日生成，使用 Julia 1.11.2；这不是文章首发日期。教程没有锁定 FiniteVolumeMethod.jl 或相关依赖的版本。

多孔介质方程描述一种非线性扩散过程。这个教程从集中在一点的初始质量出发，先估计有限时间内解可能占据的空间范围，再把无限平面上的问题近似为有限方形网格上的数值问题。第二部分加入线性增长源项，展示源项如何改变计算区域，以及如何在有限体积问题中分别传入扩散和源项参数。

用传播支撑估算计算域

无源问题写为

[
\frac{\partial u}{\partial t}=D,\nabla\cdot\left(u^\nabla u\right),
\qquad u(\mathbf,0)=M,\delta(\mathbf),
]

其中扩散函数是 (D u^)，(M) 表示初始总质量。二维 Dirac 脉冲无法直接作为网格上的普通函数取值，因此教程用窄高斯近似它：

[
g(x,y)=\frac{1}{\varepsilon^{2\pi}

\exp!\left(-\frac{x}2+y^{2}{\varepsilon}2}\right),
\qquad \varepsilon=0.1,
]

并令初始场为 (u(x,y,0)=M g(x,y))。

对这个方程，教程引用了一个有限传播支撑结果：在时刻 (t)，当

[
x^2+y2\ge R_{m,M}(Dt)^{{1/m},\qquad

R_{m,M}=\frac{4m}\left(\frac{4\pi}\right)}{(m-1)/m},
]

解为零。于是，求解到终止时刻 (T) 时，可用方形区域 (\Omega=[-L,L]^2) 代替整个平面，其中

[
L=\sqrt{R_{m,M}},(DT)^{1/(2m)}.
]

教程在该方形边界上施加零值 Dirichlet 条件。这样，区域大小来自模型参数与终止时间，而不是随意指定一个网格范围。

无源方程：从参数到有限体积求解

教程的示例参数为 (m=2)、(M=0.37)、(D=2.53)、(T=12)，并取 (\varepsilon=0.1)。先根据解析支撑计算半边长，再在正方形上生成 125×125 的三角剖分网格：

using DelaunayTriangulation, FiniteVolumeMethod

# Step 0: Define all the parameters
m = 2
M = 0.37
D = 2.53
final_time = 12.0
ε = 0.1

# Step 1: Define the mesh
RmM = 4m / (m - 1) * (M / (4π))^((m - 1) / m)
L = sqrt(RmM) * (D * final_time)^(1 / (2m))
tri = triangulate_rectangle(-L, L, -L, L, 125, 125, single_boundary=true)
mesh = FVMGeometry(tri)

原页面展示的网格摘要为 15,625 个控制体、30,752 个三角形和 46,376 条边。随后设置零值 Dirichlet 边界，把高斯初值采样到三角剖分的点上，并定义依赖状态值 (u) 的扩散函数：

# Step 2: Define the boundary conditions
BCs = BoundaryConditions(mesh, (x, y, t, u, p) -> zero(u), Dirichlet)

# Step 3: Define the actual PDE
f = (x, y) -> M * 1 / (ε^2 * π) * exp(-1 / (ε^2) * (x^2 + y^2))
diffusion_function = (x, y, t, u, p) -> p[1] * u^(p[2] - 1)
diffusion_parameters = (D, m)
initial_condition = [f(x, y) for (x, y) in DelaunayTriangulation.each_point(tri)]
prob = FVMProblem(mesh, BCs;
    diffusion_function,
    diffusion_parameters,
    initial_condition,
    final_time)

用 TRBDF2 时间积分器和 KLU 线性求解器计算，并每隔 3 个时间单位保存一次结果：

using LinearSolve, OrdinaryDiffEq
sol = solve(prob, TRBDF2(linsolve=KLUFactorization()); saveat=3.0)

页面展示的输出列出保存时刻 0.0, 3.0, 6.0, 9.0, 12.0，并把返回码显示为 Success。教程随后使用 CairoMakie 绘制时刻 0、6、12 的等值填色图：

using CairoMakie
fig = Figure(fontsize=38)
for (i, j) in zip(1:3, (1, 3, 5))
    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.005:0.05,
        colormap=:matter, extendhigh=:auto)
    tightlimits!(ax)
end
resize_to_layout!(fig)
fig

以上返回码和网格统计是原文所展示的示例输出，本稿没有运行代码，也没有独立验证误差、质量守恒或当前依赖版本兼容性。原页附有绘图示例；本稿只保留文字和绘图代码，没有复制图像。

加入线性增长源项

第二部分在扩散方程右侧加入 (\lambda u)，其中 (\lambda>0)：

[
\frac{\partial u}{\partial t}
=D,\nabla\cdot\left(u^\nabla u\right)+\lambda u.
]

初始条件仍为 (M\delta(\mathbf))，但计算域半边长按增长后的有效时间尺度调整。教程定义

[
\tau(T)=\frac{\lambda(m-1)}
\left[\exp!\left(\lambda(m-1)T\right)-1\right],
\qquad
L=\sqrt{R_{m,M}},\tau(T)^{1/(2m)}.
]

例子采用 (m=3.4)、(M=2.3)、(D=0.581)、(\lambda=0.2)、(T=10) 和 (\varepsilon=0.1)。它仍使用 125×125 网格与零值 Dirichlet 边界。代码用命名元组传递扩散参数，并把源项参数单独传入 FVMProblem。扩散函数对状态取绝对值，以便处理非整数幂：

# Step 0: Define all the parameters
m = 3.4
M = 2.3
D = 0.581
λ = 0.2
final_time = 10.0
ε = 0.1

# Step 1: Define the mesh
RmM = 4m / (m - 1) * (M / (4π))^((m - 1) / m)
L = sqrt(RmM) * (D / (λ * (m - 1)) * (exp(λ * (m - 1) * final_time) - 1))^(1 / (2m))
tri = triangulate_rectangle(-L, L, -L, L, 125, 125, single_boundary=true)
mesh = FVMGeometry(tri)

# Step 2: Define the boundary conditions
bc = (x, y, t, u, p) -> zero(u)
BCs = BoundaryConditions(mesh, bc, Dirichlet)

# Step 3: Define the actual PDE
f = (x, y) -> M * 1 / (ε^2 * π) * exp(-1 / (ε^2) * (x^2 + y^2))
diffusion_function = (x, y, t, u, p) -> p.D * abs(u)^(p.m - 1)
source_function = (x, y, t, u, λ) -> λ * u
diffusion_parameters = (D=D, m=m)
source_parameters = λ
initial_condition = [f(x, y) for (x, y) in DelaunayTriangulation.each_point(tri)]
prob = FVMProblem(mesh, BCs;
    diffusion_function,
    diffusion_parameters,
    source_function,
    source_parameters,
    initial_condition,
    final_time)

# Step 4: Solve
sol = solve(prob, TRBDF2(linsolve=KLUFactorization()); saveat=2.5)

页面输出把返回码显示为 Success，保存时刻为 0.0, 2.5, 5.0, 7.5, 10.0。绘图步骤选择 0、5、10 三个时间点：

fig = Figure(fontsize=38)
for (i, j) in zip(1:3, (1, 3, 5))
    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.05:1,
        extendlow=:auto, colormap=:matter, extendhigh=:auto)
    tightlimits!(ax)
end
resize_to_layout!(fig)
fig

两个例子共同展示了一个建模顺序：先根据方程的传播支撑估计有限域，再将初始脉冲平滑成网格可采样的窄高斯，随后依次设置网格、边界、扩散与源项、求解器和可视化。若更换参数或模型项，应重新计算域大小；示例的边长和网格分辨率不是普遍适用的固定值。

来源： Porous-Medium Equation — FiniteVolumeMethod.jl
作者： FiniteVolumeMethod.jl 文档维护者（页面未署个人名）
构建信息： Documenter.jl 1.8.0，2025-01-01；Julia 1.11.2。该信息不是文章发表日期，也不表示教程锁定了当前依赖版本。
转载许可： 原页面未列出明确的文章转载许可说明。

文章版权归作者所有，未经允许请勿转载。

THE END