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 = 0.0; # Free stream velocity
U∞ = (U,0.0);
Set up points on the body
In [4]:
n = 51;
a = 0.5; b = 0.1;
body = Bodies.Plate(1.0,n)
Out[4]:
In [6]:
plot(body)
Out[6]:
Set the motion to oscillatory pitch-heave kinematics
In [7]:
a = 0.25 # location of pitch axis, a = 0.5 is leading edge
ϕ = -π/2 # phase lag of pitch to heave
A = 0.25 # amplitude/chord
fstar = 1/π # fc/U
α₀ = 0 # mean angle of attack
Δα = 0.0 #10π/180 # amplitude of pitching
U₀ = 0.0 # translational motion (set to zero in place of free stream)
K = π*fstar # reduced frequency, K = πfc/U
oscil = RigidBodyMotions.PitchHeave(U₀,a,K,ϕ,α₀,Δα,A);
motion = RigidBodyMotion(oscil)
Out[7]:
In [60]:
motion(0.0)
Out[60]:
Transform the body with a specified initial position and orientation.
In [47]:
cent = 0.0 + 0.0im
α = 0.0
T = RigidTransform(cent,α)
T(body) # transform the body to the current configuration
Out[47]:
In [49]:
plot(body)
Out[49]:
In [11]:
xlim = (-1.0,1.0)
ylim = (-1.0,1.0)
Out[11]:
Now set up the coordinate data for operator construction
In [50]:
X = VectorData(body.x,body.y);
X̃ = VectorData(body.x̃,body.ỹ);
Set the domain size and time step size
In [53]:
Δx = 0.02;
Δt = min(0.5*Δx,0.5*Δx^2*Re);
In [63]:
sys = NavierStokes(Re,Δx,xlim,ylim,Δt,U∞ = U∞, X̃ = X̃, isstore = true, isstatic = false)
Out[63]:
Set up solution data
In [64]:
w₀ = Nodes(Dual,size(sys));
u = (w₀,[real(cent),imag(cent),α])
xg, yg = coordinates(w₀,dx=Δx,I0=Systems.origin(sys))
f = (VectorData(X̃),Vector{Float64}());
In [65]:
using LinearAlgebra
plans = ((t,u) -> Fields.plan_intfact(t,u,sys),(t,u)-> I)
Out[65]:
In [80]:
function TimeMarching.plan_constraints(u::Tuple{Nodes{Dual,NX,NY},Vector{Float64}},t,sys::NavierStokes{NX,NY,N,false}) where {NX,NY,N}
# for now, just assume that there is only one body. will fix later.
xc, yc, α = u[2]
T = Bodies.RigidTransform((xc,yc),α)
# should be able to save some time and memory allocation here...
x, y = T(sys.X̃.u,sys.X̃.v)
X = VectorData(x,y)
regop = Regularize(X,cellsize(sys);issymmetric=true,I0=origin(sys))
if sys._isstore
Hmat, Emat = RegularizationMatrix(regop,VectorData{N}(),Edges{Primal,NX,NY}())
sys.Hmat = Hmat
sys.Emat = Emat
return (f->TimeMarching.B₁ᵀ(f,sys), f->zeros(Float64,size(u[2]))),
(w->TimeMarching.B₂(w,sys), u->Vector{Float64}())
else
return (f->TimeMarching.B₁ᵀ(f,regop,sys), f->zeros(Float64,size(u[2]))),
(w->TimeMarching.B₂(w,regop,sys), u->Vector{Float64}())
end
end
In [82]:
ifherk = IFHERK(u,f,sys.Δt,
plans,
(u,t) -> TimeMarching.plan_constraints(u,t,sys),
((u,t) -> Systems.r₁(u,t,sys,motion),
(u,t) -> Systems.r₂(u,t,sys,motion)),tol=1e-3,rk=TimeMarching.RK31,isstored=true,isstaticconstraints=false)
Out[82]:
In [83]:
t = 0.0
w₀ .= 0.0
u = (deepcopy(w₀),[real(cent),imag(cent),α])
f = (VectorData(X̃),Vector{Float64}());
tf = 0.1;
T = Δt:Δt:tf;
fx = Float64[];
fy = Float64[];
thist = [];
uhist = [];
tsample = 0.1;
In [86]:
@time for ti in T
global t, u, f = ifherk(t,u)
# save data for later use
push!(thist,t)
push!(fx,sum(f[1].u)*Δx^2)
push!(fy,sum(f[1].v)*Δx^2)
(isapprox(mod(t,tsample),0,atol=1e-12) || isapprox(mod(t,tsample),tsample,atol=1e-12)) ? push!(uhist,u) : nothing
end
println("solution completed through time t = ",t)
In [87]:
plot(xg,yg,uhist[1][1],levels=range(-1,1,length=30),color = :RdBu)
Tr = RigidTransform(uhist[1][2])
Tr(body) # transform the body to the current configuration
plot!(body,linecolor=:black)
Out[87]:
In [88]:
px = plot(thist,2*fx,ylim=(-3,3),xlabel="Convective time",ylabel="\$C_D\$",legend=false)
py = plot(thist,2*fy,ylim=(-6,6),xlabel="Convective time",ylabel="\$C_L\$",legend=false)
plot(px,py)
Out[88]:
In [90]:
ψ = deepcopy(u[1])
ψ .= -(sys.L\u[1])*Δx .+ repeat(collect(yg)',size(sys,1),1)
plot(xg,yg,ψ,trim=1,clim=(-0.1,0.1))
plot!(body,fillcolor=:black,linecolor=:black)
Out[90]:
In [ ]: