这篇文章将演示如何使用 Julia 和 WaterLily.jl 设置并模拟一个游泳的角鲨模型。这篇文章改编自一个 Pluto 笔记本，你可以从 这里 获取。
using WaterLily, StaticArrays, PlutoUI, Interpolations, Plots, Images我们将使用一个简单的鲨鱼模型，该模型基于 Lighthill 关于 细长鱼类游泳的开创性论文。它专注于鱼的“脊骨”；将形状理想化为中心两侧的厚度分布，并将运动理想化为横向（“左右”）行进波。令人惊讶的是，这种简单的方法提供了对广泛的游泳动物的洞察力，如下图所示。 图片来自：Gazzola 等人，自然杂志 2014 年
我们将使用插值来定义厚度和运动分布，以便通过一些点拟合样条曲线。
fit = y -> scale(
interpolate(y, BSpline(Quadratic(Line(OnGrid())))),
range(0,1,length=length(y))
)这是一张我们将模拟的角鲨的图片。图片来自：鲨鱼信托基金。
url2 = "https://pterosaurheresies.files.wordpress.com/2020/01/squalus-acanthias-invivo588.jpg"
filename2 = download(url2)
dogfish = load(filename2)
底部视图显示了我们感兴趣的轮廓，沿着长度添加一些点定义了厚度分布函数 thk。
plot(dogfish)
nose, len = (30,224), 500
width = [0.02, 0.07, 0.06, 0.048, 0.03, 0.019, 0.01]
scatter!(
nose[1] .+ len .* range(0, 1, length=length(width)),
nose[2] .- len .* width, color=:blue, legend=false)
thk = fit(width)
x = 0:0.01:1
plot!(
nose[1] .+ len .* x,
[nose[2] .- len .* thk.(x), nose[2] .+ len .* thk.(x)],
color=:blue)
观察 游泳的角鲨的视频，我们可以看到一些一般特征
身体前半部分的运动幅度很小（大约是尾巴的 20%）。这为行进波设定了幅度包络。
envelope = [0.2,0.21,0.23,0.4,0.88,1.0]
amp = fit(envelope)行进波的波长略大于身体长度
在下面的代码中，您可以更改 的值来控制行进波的波长，观察它对运动周期内脊柱的影响。
λ = 1.4
scatter(0:0.2:1, envelope)
colors = palette(:cyclic_wrwbw_40_90_c42_n256)
for t in 1/12:1/12:1
plot!(x, amp.(x) .* sin.(2π/λ * x .- 2π*t),
color=colors[floor(Int,t*256)])
end
plot!(ylim=(-1.4,1.4), legend=false)
现在厚度和运动已经定义了，但我们如何将它们应用于流体模拟？WaterLily 使用 浸入边界法 和 自动微分 将物体嵌入流体中。结果是，我们不需要进行任何网格划分；我们只需要一个到表面的有符号距离函数 (SDF)。
让我们从定义到“脊柱”的 SDF 开始，这是一个从 的线段。 观看来自 Inigo Quilez 的这个很棒的视频，了解这种 sdf 的推导。。下面的图显示了 sdf 和零轮廓，它仅仅是一条线段。
对 SDF 进行简单的调整可以让我们更好地控制形状和位置。通过将 y 偏移量更改为 y = y-shift，我们可以横向移动物体。通过从距离中减去厚度，例如 sdf = sdf-thickness，我们可以给线条一些宽度。这就是模拟鲨鱼所需的一切。
shift = 0.5
T = 0.5
function segment_sdf(x, y)
s = clamp(x, 0, 1) # distance along the segment
y = y - shift # shift laterally
sdf = √sum(abs2, (x-s, y)) # line segment SDF
return sdf - T * thk(s) # subtract thickness
end
grid = -1:0.05:2
contourf(grid, grid, segment_sdf, clim=(-1,2), linewidth=0)
contour!(grid, grid, segment_sdf, levels=[0], color=:black) # zero contour
在测试了基本的 SDF 之后，我们就可以使用下面定义的 fish 函数来设置 WaterLily 模拟。
thk 函数用于创建 sdf，amp 函数用于创建行进波 map。
传递到 fish 的唯一数值参数是鱼的长度 L，以计算单元为单位。这会设置模拟的分辨率和流体数组的大小。
function fish(thk, amp, k=5.3; L=2^6, A=0.1, St=0.3, Re=1e4)
# fraction along fish length
s(x) = clamp(x[1]/L, 0, 1)
# fish geometry: thickened line SDF
sdf(x,t) = √sum(abs2, x - L * SVector(s(x), 0.)) - L * thk(s(x))
# fish motion: travelling wave
U = 1
ω = 2π * St * U/(2A * L)
function map(x, t)
xc = x .- L # shift origin
return xc - SVector(0., A * L * amp(s(xc)) * sin(k*s(xc)-ω*t))
end
# make the fish simulation
return Simulation((4L+2,2L+2), [U,0.], L;
ν=U*L/Re, body=AutoBody(sdf,map))
end
# Create the swimming shark
L,A,St = 3*2^5,0.1,0.3
swimmer = fish(thk, amp; L, A, St);
# Save a time span for one swimming cycle
period = 2A/St
cycle = range(0, 23/24*period, length=24)我们可以通过绘制浸入边界函数 μ₀ 来测试几何形状；它在流体中等于 1，在物体中等于 0。
@gif for t ∈ cycle
measure!(swimmer, t*swimmer.L/swimmer.U)
contour(swimmer.flow.μ₀[:,:,1]',
aspect_ratio=:equal, legend=false, border=:none)
end
运动的动画看起来很棒，所以我们准备运行流体模拟器了！
sim_step!(sim, t, remeasure=true) 函数运行模拟器到时间 t，在每次步长时重新测量物体的位置。（默认情况下，remeasure=false，因为这需要一些额外的计算时间，并且对于静态几何形状不需要。）
# run the simulation a few cycles (this takes few seconds)
sim_step!(swimmer, 10, remeasure=true)
sim_time(swimmer)模拟现在已经向前运行了一段时间，但默认情况下没有可视化或测量。
为了观察发生的事情，让我们制作一个涡量 ω=curl(u) 的 GIF 来可视化鲨鱼尾迹中的涡流。这需要模拟一个运动周期，并在模拟中计算所有 @inside 点的 curl。
# plot the vorcity ω=curl(u) scaled by the body length L and flow speed U
function plot_vorticity(sim)
@inside sim.flow.σ[I] = WaterLily.curl(3, I, sim.flow.u) * sim.L / sim.U
contourf(sim.flow.σ',
color=palette(:BuGn), clims=(-10, 10), linewidth=0,
aspect_ratio=:equal, legend=false, border=:none)
end
# make a gif over a swimming cycle
@gif for t ∈ sim_time(swimmer) .+ cycle
sim_step!(swimmer, t, remeasure=true)
plot_vorticity(swimmer)
end
这很漂亮（毕竟 CFD 代表的是彩色流体动力学），但也告诉我们关于流动的重要信息。请注意，除了尾巴之外，物体上没有任何其他地方出现涡流！这表明效率很高，因为能量只用于产生这些尾部涡流。
我们也可以从模拟中获取一些定量测量。∮nds 函数在物体表面上进行积分。通过传入压力 p，我们可以测量鲨鱼产生的推力和侧向力！
function get_force(sim, t)
sim_step!(sim, t, remeasure=true)
return WaterLily.∮nds(sim.flow.p, sim.body, t*sim.L/sim.U) ./ (0.5*sim.L*sim.U^2)
end
forces = [get_force(swimmer, t) for t ∈ sim_time(swimmer) .+ cycle]
scatter(cycle./period, [first.(forces), last.(forces)],
labels=permutedims(["thrust", "side"]),
xlabel="scaled time",
ylabel="scaled force")
我们可以从这张简单的图中了解到很多东西。例如，侧向力的频率与游泳运动本身的频率相同，而推力的频率是其两倍，每次尾巴经过中心线时都会出现一个峰值。
这个简单的模型是一个很好的起点，它为改进鲨鱼模拟和提出研究问题开辟了许多途径
在自由游泳物体中，瞬时净力应该为零！我们可以向我们的模型中添加反应运动来实现这一点。如果我们这样做，模型鲨鱼会直线游泳吗？还是需要控制回路？
真实的鲨鱼是 3D 的（哎呀！）。虽然我们可以轻松地使用 2D 样条曲线扩展这种方法，但模拟所需的时间会长得多。有没有办法使用 GPU 来加速模拟，而无需完全改变求解器？
如果我们要制造一个受这种鲨鱼启发的仿生机器人，我们将对形状、运动和可用动力有约束。我们可以使用此框架来帮助在约束条件下优化我们的机器人吗？
在下面，您可以找到在本笔记本中使用的所有软件包的链接。祝您模拟愉快！
