In [1]:
using ViscousFlow
In [2]:
using Plots
pyplot()
clibrary(:colorbrewer)
default(grid = false)
Set the flow parameters
In [3]:
Re = 200; # Reynolds number
U = 1.0; # Free stream velocity
U∞ = (U,0.0);
Set up points on the body.
In [4]:
n = 100;
body = Bodies.Ellipse(0.5,n)
Out[4]:
Transform the body to lie at a specific position. Here, we just put it at the origin.
In [5]:
cent = (0.0,0.0)
α = 0.0
T = RigidTransform(cent,α)
T(body) # transform the body to the current configuration
Out[5]:
Set up the domain and plot the body inside it.
In [6]:
xlim = (-1.0,3.0)
ylim = (-1.0,1.0);
In [7]:
plot(body,xlim=xlim,ylim=ylim)
Out[7]:
Now set up the coordinate data for operator construction
In [8]:
X = VectorData(body.x,body.y);
#X̃ = VectorData(body.x̃,body.ỹ);
Set the grid cell size and time step size.
In [9]:
Δx = 0.02
cfl = 0.5
Δt = min(cfl*Δx,0.5*Δx^2*Re)
Out[9]:
Set up the state vector and constraint force vector for a static body
In [10]:
sys = Systems.NavierStokes(Re,Δx,xlim,ylim,Δt,U∞ = U∞, X̃ = X, isstore = true)
Out[10]:
In [11]:
w₀ = Nodes(Dual,size(sys));
f = VectorData(X);
In [12]:
wf = Systems.PointForce(w₀,(1.5,0.0),10.0,4.0,1.0,sys)
Out[12]:
In [13]:
xg, yg = coordinates(w₀,dx=Δx,I0=Systems.origin(sys))
Out[13]:
Set up the integrator here
In [14]:
plan_intfact(t,u) = Systems.plan_intfact(t,u,sys)
plan_constraints(u,t) = TimeMarching.plan_constraints(u,t,sys)
r₁(u,t) = TimeMarching.r₁(u,t,sys) + wf(t)
r₂(u,t) = TimeMarching.r₂(u,t,sys)
@time ifherk = IFHERK(w₀,f,sys.Δt,plan_intfact,plan_constraints,(r₁,r₂),rk=TimeMarching.RK31,isstored=true)
Out[14]:
Initialize the state vector and the history vectors
In [15]:
t = 0.0
u = deepcopy(w₀)
fx = Float64[];
fy = Float64[];
thist = Float64[];
uhist = [];
tsample = 0.2; # rate at which to store field data
Set the time range to integrate over.
In [16]:
tf = 30.0;
T = Δt:Δt:tf;
In [17]:
@time for ti in T
global t, u, f = ifherk(t,u)
push!(thist,t)
push!(fx,sum(f.u)*Δx^2)
push!(fy,sum(f.v)*Δx^2)
(isapprox(mod(t,tsample),0,atol=1e-8) || isapprox(mod(t,tsample),tsample,atol=1e-8)) ? push!(uhist,u) : nothing
end
println("solution completed through time t = ",t)
In [18]:
thist[end]
Out[18]:
Basic plot
In [19]:
plot(xg,yg,uhist[end],levels=range(-0.1,stop=0.1,length=30), color = :RdBu,xlim=(-1+Δx,3-Δx),ylim=(-1+Δx,1-Δx),width=1)
plot!(body)
Out[19]:
Make a movie
In [20]:
anim = @animate for i = 1:length(uhist)
plot(xg,yg,uhist[i],levels=range(-0.1,stop=0.1,length=30),color = :RdBu,xlim=(-1+Δx,3-Δx),ylim=(-1+Δx,1-Δx))
plot!(body)
end
gif(anim,"cylinderRe200.gif")
Out[20]:
In [21]:
plt = plot(layout = (2,1), size = (600, 400))
plot!(plt[1],thist,2*fy,xlim=(0,30),ylim=(-2,2),xlabel="Convective time",ylabel="\$C_L\$",legend=false)
plot!(plt[2],thist,2*fx,xlim=(0,30),ylim=(0,4),xlabel="Convective time",ylabel="\$C_D\$",legend=false)
plt
Out[21]:
In [22]:
using Statistics
In [23]:
thist[3000]
Out[23]:
In [30]:
Statistics.mean(2*fx)
Out[30]:
In [32]:
using Pkg
In [33]:
Pkg.add("FFTW")
In [34]:
using FFTW
In [35]:
fft(fx)
Out[35]:
In [36]:
nt = length(thist)
nfreq = (nt-1)÷2+1
fsamp = 1/Δt
freq = 0.5*fsamp*range(0,stop=1.0,length=nfreq)
fhat = fft(fy[1:nt-1])/(nt-1)
plot(freq,4*abs.(fhat[1:nfreq]),xlim=(0,1))
Out[36]:
In [39]:
xg,yg = coordinates(u,dx=Δx,I0=Systems.origin(sys))
ψ = deepcopy(uhist[end])
ψ .= -(sys.L\uhist[end])*Δx .+ repeat(collect(yg)',size(sys,1))
plot(xg,yg,ψ,trim=1,clim=(-0.1,0.1))
plot!(body,fillcolor=:black,linecolor=:black)
Out[39]:
In [ ]: