二维鱼游CFD代码

最近学了会Julia，参考了原作者的shark，做一下基于airfoils 2D的鱼游，暂时没想好有什么需要深入研究的，代码公开如下：

鱼身是naca0016，然后一些参数可以参考我以前发的论文。

using WaterLily, StaticArrays, Interpolations, Plots, Images

fit = y -> scale(interpolate(y, BSpline(Quadratic(Line(OnGrid())))),range(0,1,length=length(y)))
width1 = [0.02, 0.07, 0.06, 0.048, 0.03, 0.019, 0.01]
function naca0016_points(n_points::Int)
    x = range(0, stop=1, length=n_points)
    yt = 5 * 0.16 * (0.2969 * sqrt.(x) .- 0.1260 * x .- 0.3516 * x.^2 .+ 0.2843 * x.^3 .- 0.1015 * x.^4)
    return yt
end
width2 = naca0016_points(100)
thk = fit(width2)
envelope1 = [0.2,0.21,0.23,0.4,0.88,1.0]
ampfun = (x; a0=0, a1=0.00294, a2=0.111) -> a0 + a1 * x + a2 * x^2
envelop2 = ampfun.(0:0.01:1; a0=0, a1=0.00294, a2=0.111)
amp = fit(envelop2)

function fish(thk, amp, k=10.472;St,ω, L=2^6,fish_shift, Re=1e4, U = 1,T=Float32,mem=Array)
  function map(x, t)
		xc = x .- fish_shift # shift origin
		return xc - SVector(0., L * amp(s(xc)) * sin(k*s(xc)-ω*t))
	end
	s(x) = clamp(x[1]/L, 0, 1)
	sdf(x,t) = √sum(abs2, x - L * SVector(s(x), 0.)) - L * thk(s(x))
	return Simulation((300,180), (U,0.), L; ν=U*L/Re, body=AutoBody(sdf,map))
end
pending_singal = (x, coeff) -> coeff * x

U = 1
L = 60 # length of the fish L
fish_shift = 90
St = 6.858*amp(1)/U
ω = 2π * St * U/(2amp(1) * L)
swimmer = fish(thk, amp;St,ω, L, U, fish_shift);

period = 2amp(1)/St # Save a time span for one swimming cycle
cycle = range(0, 23/24*period, length=24)
@gif for t ∈ cycle
	measure!(swimmer, t*swimmer.L/swimmer.U)
	contour(swimmer.flow.μ₀[:,:,1]',
			aspect_ratio=:equal, legend=false, border=:none)
end

结果可以看到鱼游动的形态

在这里插入图片描述

下面是仿真

stop = 5
@time @gif for tᵢ in range(0.,stop;step=0.05)
    if tᵢ % 0.2 == 0
      println("Simulation tU/L = ",round(tᵢ,digits=4))
    end
    sim_step!(swimmer,tᵢ,remeasure=true)
end

保存gif

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
cycle = range(0, 23/24*period, length=24)
# make a gif over a swimming cycle
@gif for t ∈ sim_time(swimmer) .+ cycle
	sim_step!(swimmer, t, remeasure=true)
	plot_vorticity(swimmer)
end