In [1]:
using ViscousFlow
In [2]:
import ViscousFlow.Fields: curl
include("./Streaming.jl")
In [3]:
using Plots
pyplot()
clibrary(:colorbrewer)
default(grid = false)
In [4]:
using Statistics
The motion of the body is described by oscillatory small-amplitude translation,
$V_b(\xi,t) = \epsilon \hat{U}(t)$
where $\epsilon \ll 1$, $\hat{U}(t)$ is a periodic velocity with unit amplitude, and $\xi \in S_b$ (i.e. it lies on the body surface). The associated position of any point on the surface of the body is then described by
$X_b(\xi,t) = \xi + \epsilon \int_0^t \hat{U}(\tau)d\tau.$
We will write the flow field in asympotic form, e.g. $v = \epsilon v_1 + \epsilon^2 v_2$, and seek to solve two asymptotic levels of equation (for vorticity)
$\dfrac{\partial w_1}{\partial t} - \dfrac{1}{Re} \nabla^2 w_1 = 0$
subject to boundary condition $v_1(\xi,t) = \hat{U}(t)$ for all $\xi \in S_b$. Note that the boundary condition is applied at the initial location of the surface, not its time-varying location.
And at second order,
$\dfrac{\partial w_2}{\partial t} - \dfrac{1}{Re} \nabla^2 w_2 = -\nabla\cdot(v_1 w_1),$
subject to boundary condition $v_2(\xi,t) = -\int_0^t \hat{U}(\tau)d\tau \cdot \nabla v_1(\xi,t)$ for all $\xi \in S_b$. This is also applied at the initial location of the surface.
Thus, to solve this problem, we will set up a state vector that contains $w_1$, $w_2$, and an unscaled 'body' state $x_c$, $y_c$. These latter components will be used to hold the components of $\Delta\hat{X} \equiv \int_0^t \hat{U}(\tau)d\tau$.
A fluid particle, initially at a location $x$, is subjected to a slightly different velocity, $U(x,t)$, over the ensuing oscillation period, since it is advected by the local velocity to sample nearby points. Its first-order velocity is simply $U_1(x,t) = v_1(x,t)$. Its second order velocity, however, is
$U_2(x,t) = v_2(x,t) + \int_0^t v_1(x,\tau)d\tau \cdot \nabla v_1(x,t)$
In particular, if the particle is initially on the surface, $x \in S_b$, then $v_1 = \hat{U}$, and the second term of this expression above cancels $v_2$, so that $U_2 = 0$. This ensures that the particle continues to oscillate with the body's surface velocity. Let us define an Eulerian displacement field,
$\Delta X(x,t) = \int_0^t v(x,\tau)d\tau, \quad \Delta X(x,0) = 0,$
or, equivalently,
$\dfrac{d\Delta X}{dt} = v(x,t), \quad \Delta X(x,0) = 0.$
Then, expanding $\Delta X$ in an asymptotic sequence, $\Delta X = \epsilon \Delta X_1 + \epsilon^2 \Delta X_2$, we have
$U_2(x,t) = v_2(x,t) + \Delta X_1(x,t) \cdot \nabla v_1(x,t)$
Let us define a drift displacement, $\Delta X_d$, by the solution of
$\dfrac{d\Delta X_d}{dt} = \Delta X_1(x,t) \cdot \nabla v_1(x,t), \quad \Delta X_d(x,0) = 0$
The motion of a fluid particle, starting at $x$, is thus given by
$X_p(x,t) = x + \epsilon \Delta X_1(x,t) + \epsilon^2 (\Delta X_2(x,t) + \Delta X_d(x,t))$
These represent the integral curves of the velocity field. If we look only at these integral curves over one period, then we seek the net displacement made by a particle over this period,
$X_p(x,T) - x = \Delta X(x,T) + \epsilon^2 \Delta X_d(x,T) = \int_0^T U(x,\tau)d\tau \equiv T \overline{U}(x).$
We only need a mean velocity field for this purpose.
In [5]:
Re = 40.0
ϵ = 0.1
Ω = 1.0 # angular frequency
Ax = 1.0 # x amplitude (before scaling by ϵ)
ϕx = 0.0 # phase lead of x motion
Ay = 0.0 # y amplitude
ϕy = 0.0 # phase lead of y motion
oscil = RigidBodyMotions.RigidBodyMotion(RigidBodyMotions.Oscillation(Ω,Ax,ϕx,Ay,ϕy))
Out[5]:
In [6]:
t = range(0.0,stop=4π,length=401)
u = map(ti -> real(oscil(ti)[2]),t) # u component of velocity of centroid
plot(t,u)
Out[6]:
In [7]:
p = Params(0.1,Re)
Out[7]:
First-order solution
In [8]:
K = 1
Ψ₁ = ComplexFunc(r -> -p.C/r + 2Y(r)/p.γ)
W₁ = D²(Ψ₁,K) # note that this is actually the negative of the vorticity. We will account for this when we evaluate it.
Ur₁, Uθ₁ = curl(Ψ₁,K)
# for verifying the solution
LW₁ = D²(W₁,K);
resid1 = ComplexFunc(r -> LW₁(r)+im*p.Re*W₁(r))
println("Maximum residual on solution = ",maximum(abs.(resid1.(range(1,5,length=10)))))
# for verifying boundary conditions
dΨ₁ = ComplexFunc(r -> derivative(Ψ₁,r))
bcresid1 = Ψ₁(1) - 1
bcresid2 = dΨ₁(1) - 1
println("BC residual on Ψ₁(1) = ",abs(bcresid1))
println("BC residual on dΨ₁(1) = ",abs(bcresid2))
In [9]:
function ω₁ex(x,y,t)
r = sqrt(x^2+y^2)
return real(-W₁(r)*y/r*exp.(-im*t))
end
function u₁ex(x,y,t)
r = sqrt(x^2+y^2)
coseval = x/r
sineval = y/r
return real.((Ur₁(r)*coseval^2-Uθ₁(r)*sineval^2)*exp.(-im*t))
end
function v₁ex(x,y,t)
r = sqrt(x^2+y^2)
coseval = x/r
sineval = y/r
return real.((Ur₁(r)+Uθ₁(r))*coseval*sineval*exp.(-im*t))
end
function ψ₁ex(x,y,t)
r = sqrt(x^2+y^2)
return real(Ψ₁(r)*y/r*exp.(-im*t))
end
Out[9]:
Analytical solution of the steady part
In [10]:
K = 2
fakefact = 1
#f₀ = ComplexFunc(r -> -0.5*p.γ²*p.Re*(0.5*(p.C*conj(X(r))-conj(p.C)*X(r))/r^2 + X(r)*conj(Z(r)) - conj(X(r))*Z(r)))
f₀ = ComplexFunc(r -> -p.γ²*p.Re*(0.5*p.C*conj(X(r))/r^2 + X(r)*conj(Z(r))))
f̃₀ = ComplexFunc(r -> f₀(r) - 0.5*p.γ²*p.Re*(-0.5*conj(Z(r))+0.5*Z(r)))
I⁻¹ = ComplexIntegral(r->f₀(r)/r,1,Inf,length=100000)
I¹ = ComplexIntegral(r->f₀(r)*r,1,Inf,length=100000)
I³ = ComplexIntegral(r->f₀(r)*r^3,1,20,length=400000)
I⁵ = ComplexIntegral(r->f₀(r)*r^5,1,20,length=400000)
Ψs₂ = ComplexFunc(r -> -r^4/48*I⁻¹(r) + r^2/16*I¹(r) + I³(r)/16 + I⁻¹(1)/16 - I¹(1)/8 - fakefact*0.25im*p.γ*Y(1) +
1/r^2*(-I⁵(r)/48-I⁻¹(1)/24+I¹(1)/16 + fakefact*0.25im*p.γ*Y(1)))
Ws₂ = D²(Ψs₂,K)
Usr₂, Usθ₂ = curl(Ψs₂,K)
# for verifying the solution
LWs₂ = D²(Ws₂,K)
resids = ComplexFunc(r -> real(LWs₂(r)-f₀(r)))
println("Maximum residual on solution = ",maximum(abs.(resids.(range(1,5,length=10)))))
# for verifying boundary conditions
dΨs₂ = ComplexFunc(r -> derivative(Ψs₂,r))
bcresids1 = Ψs₂(1)
bcresids2 = real(dΨs₂(1) - 0.25im*W₁(1))
println("BC residual on Ψs₂(1) = ",abs(bcresids1))
println("BC residual on dΨs₂(1) = ",abs(bcresids2))
Analytical solution of unsteady part
In [11]:
K = 2
fakefact = 1
g₀ = ComplexFunc(r -> 0.5*p.γ²*p.Re*p.C*X(r)/r^2)
g̃₀ = ComplexFunc(r -> g₀(r) - 0.5*p.γ²*p.Re*Z(r))
Kλ = ComplexFunc(r -> H11(1)*H22(r) - H12(1)*H21(r))
IKgr = ComplexIntegral(r -> r*Kλ(r)*g₀(r),1,20,length=400000)
IH21gr = ComplexIntegral(r -> r*H21(r)*g₀(r),1,Inf,length=100000)
Igr⁻¹ = ComplexIntegral(r -> g₀(r)/r,1,Inf,length=100000)
Igr³ = ComplexIntegral(r -> g₀(r)*r^3,1,20,length=400000)
Ig¹ = ComplexFunc(r -> 0.25im*π/(p.λ²*H11(1))*IKgr(r)*H21(r))
Ig² = ComplexFunc(r -> 0.25im*π/(p.λ²*H11(1))*IH21gr(r)*Kλ(r))
Ig³ = ComplexFunc(r -> 1/(p.λ²*p.λ*H11(1))*((H21(r)-H21(1)/r^2)*Igr⁻¹(1)+IH21gr(1)/r^2))
Ig⁴ = ComplexFunc(r -> -0.25/p.λ²*(Igr⁻¹(r)*r^2-Igr⁻¹(1)/r^2+Igr³(r)/r^2))
Ψ₂ = ComplexFunc(r -> Ig¹(r) + Ig²(r) + Ig³(r) + Ig⁴(r) + fakefact*0.5im/sqrt(2)*Y(1)/H11(1)*(H21(r)-H21(1)/r^2))
Ψ̃₂ = ComplexFunc(r -> Ψ₂(r)+ 0.5im*(-p.C/r^2 + Z(r))) # cylinder-fixed reference frame... not used
W₂ = D²(Ψ₂,K)
Ur₂, Uθ₂ = curl(Ψ₂,K);
# for verifying the solution
LW₂ = D²(W₂,K);
resid = ComplexFunc(r -> LW₂(r)+2im*p.Re*W₂(r)-g₀(r))
println("Maximum residual on solution = ",maximum(abs.(resid.(range(1,5,length=10)))))
# for verifying boundary conditions
dΨ₂ = ComplexFunc(r -> derivative(Ψ₂,r))
bcresid1 = Ψ₂(1)
bcresid2 = dΨ₂(1) - 0.5im*p.γ*Y(1)
println("BC residual on Ψ₂(1) = ",abs(bcresid1))
println("BC residual on dΨ₂(1) = ",abs(bcresid2))
In [12]:
function ω₂ex(x,y,t)
r = sqrt(x^2+y^2)
sin2eval = 2*x*y/r^2
return real(-Ws₂(r) .- W₂(r)*exp.(-2im*t))*sin2eval
end
function u₂ex(x,y,t)
r = sqrt(x^2+y^2)
coseval = x/r
sineval = y/r
cos2eval = coseval^2-sineval^2
sin2eval = 2*coseval*sineval
ur = real.(Usr₂(r) .+ Ur₂(r)*exp.(-2im*t))*cos2eval
uθ = real.(Usθ₂(r) .+ Uθ₂(r)*exp.(-2im*t))*sin2eval
return ur*coseval .- uθ*sineval
end
function v₂ex(x,y,t)
r = sqrt(x^2+y^2)
coseval = x/r
sineval = y/r
cos2eval = coseval^2-sineval^2
sin2eval = 2*coseval*sineval
ur = real.(Usr₂(r) .+ Ur₂(r)*exp.(-2im*t))*cos2eval
uθ = real.(Usθ₂(r) .+ Uθ₂(r)*exp.(-2im*t))*sin2eval
return ur*sineval .+ uθ*coseval
end
function ψ₂ex(x,y,t)
r = sqrt(x^2+y^2)
coseval = x/r
sineval = y/r
sin2eval = 2*coseval*sineval
return real(Ψs₂(r) .+ Ψ₂(r)*exp.(-2im*t))*sin2eval
end
function ψ̄₂ex(x,y)
r = sqrt(x^2+y^2)
coseval = x/r
sineval = y/r
sin2eval = 2*coseval*sineval
return real(Ψs₂(r))*sin2eval
end
Out[12]:
In [17]:
n = 150;
body = Ellipse(1.0,n)
Out[17]:
Transform the body with a specified initial position and orientation.
In [14]:
cent = (0.0,0.0)
α = 0.0
T = RigidTransform(cent,α)
T(body) # transform the body to the current configuration
Out[14]:
Now set up the coordinate data for operator construction
In [15]:
coords = VectorData(body.x,body.y);
In [18]:
xlim = (-4.0,4.0)
ylim = (-4.0,4.0)
Out[18]:
In [19]:
plot(body,xlim=xlim,ylim=ylim)
Out[19]:
Set the cell size and time step size
In [32]:
Δx = 0.0202
Δt = min(π*0.5*Δx,π*0.5*Δx^2*Re)
Out[32]:
The discretized equations include the constraint forces on the surface, used to enforce the boundary conditions. The first level can be written as
$\dfrac{dw_1}{dt} - \dfrac{1}{Re}Lw_1 + C^T E^T f_1 = 0,\quad -E C L^{-1}w_1 = \hat{U}, \quad w_1(0) = 0$
where $E$ is the interpolation operator, $E^T$ is the regularization operator, $C^T$ is the rot (curl) operator, and $f_1$ represents the vector of discrete surface forces.
The second asymptotic level is written as
$\dfrac{dw_2}{dt} - \dfrac{1}{Re}Lw_2 + C^T E^T f_2 = -N(w_1),\quad -E C L^{-1}w_2 = \Delta\hat{X}\cdot E G C L^{-1}w_1, \quad w_2(0) = 0$
where $N$ represents the non-linear convective term. We must also advance the state $\Delta\hat{X}$:
$\dfrac{d\Delta\hat{X}}{dt} = \hat{U}, \quad \Delta\hat{X}(0) = 0$
To account for the fluid particle motion, we will also integrate the Eulerian displacement fields $\Delta X_1$ and $\Delta X_2$:
$\dfrac{d\Delta X_1}{dt} = v_1, \quad \Delta X_1(x,0) = 0$
and
$\dfrac{d\Delta X_2}{dt} = v_2, \quad \Delta X_2(x,0) = 0$
The construction of the NavierStokes system creates and stores the regularization and interpolation matrices. These will not change with time, since the boundary conditions are applied at the initial location of the surface points.
In [33]:
ddftype = Fields.Goza
sys = NavierStokes(Re,Δx,xlim,ylim,Δt, X̃ = coords, isstore = true, isasymptotic = true, isfilter = false, ddftype = ddftype)
Out[33]:
Set up the state vector and constraint force vector for a static body
In [34]:
w₀ = Nodes(Fields.Dual,size(sys))
ΔX = Edges(Primal,w₀)
f1 = VectorData(coords);
xg, yg = coordinates(w₀,sys.grid)
Out[34]:
Make the tuple data structures. The state tuple holds the first and second asymptotic level vorticity, the Eulerian displacement field for each level, and the two components of the unscaled centroid displacement of the body. The last three parts have no constraints, so we set the force to empty vectors.
In [35]:
u = (zero(w₀),zero(w₀),zero(ΔX),zero(ΔX),[0.0,0.0])
f = (zero(f1),zero(f1),Vector{Float64}(),Vector{Float64}(),Vector{Float64}())
TU = typeof(u)
TF = typeof(f)
Out[35]:
The right-hand side of the first-order equation is 0. The right-hand side of the second-order equation is the negative of the non-linear convective term, based on the first-order solution; for this, we use the predefined r₁. The right-hand side of the Eulerian displacement field equations are just the corresponding velocities at those levels. The right-hand side of the body update equation is the unscaled velocity.
In [36]:
function TimeMarching.r₁(u::TU,t::Float64)
_,ċ,_,_,α̇,_ = oscil(t)
return zero(u[1]), TimeMarching.r₁(u[1],t,sys), lmul!(-1,curl(sys.L\u[1])), lmul!(-1,curl(sys.L\u[2])), [real(ċ),imag(ċ)]
end
The right-hand side of the first constraint equation is the unscaled rigid-body velocity, evaluated at the surface points. The right-hand side of the second constraint equation is $-\hat{X}\cdot\nabla v_1$. The Eulerian displacement and the body update equations have no constraint, so these are set to an empty vector.
In [37]:
gradopx = Regularize(sys.X̃,cellsize(sys);I0=origin(sys),issymmetric=true,ddftype=ddftype,graddir=1)
gradopy = Regularize(sys.X̃,cellsize(sys);I0=origin(sys),issymmetric=true,ddftype=ddftype,graddir=2)
dEx = InterpolationMatrix(gradopx,sys.Fq,sys.Vb)
dEy = InterpolationMatrix(gradopy,sys.Fq,sys.Vb)
nothing
In [38]:
function TimeMarching.r₂(u::TU,t::Float64)
fact = 2 # not sure how to explain this factor yet.
_,ċ,_,_,α̇,_ = oscil(t)
U = (real(ċ),imag(ċ))
# -ΔX̂⋅∇v₁
Δx⁻¹ = 1/cellsize(sys)
sys.Fq .= curl(sys.L\u[1]) # -v₁
sys.Vb .= dEx*sys.Fq # dv₁/dx*Δx
sys.Vb.u .*= -fact*Δx⁻¹*u[5][1]
sys.Vb.v .*= -fact*Δx⁻¹*u[5][1]
Vb = deepcopy(sys.Vb) # -X⋅dv₁/dx
sys.Vb .= dEy*sys.Fq # dv₁/dy*Δx
sys.Vb.u .*= -fact*Δx⁻¹*u[5][2]
sys.Vb.v .*= -fact*Δx⁻¹*u[5][2]
Vb .+= sys.Vb # -X⋅dv₁/dx - Y.⋅dv₁/dy
return U + α̇ × (sys.X̃ - body.cent), Vb, Vector{Float64}(), Vector{Float64}(), Vector{Float64}()
end
The integrating factor for the first two equations is simply the one associated with the usual viscous term. The last three equations have no term that needs an integrating factor, so we set their integrating factor operators to the identity, I.
In [39]:
using LinearAlgebra
plans = ((t,u) -> Fields.plan_intfact(t,u,sys),(t,u) -> Fields.plan_intfact(t,u,sys),(t,u) -> Identity(),(t,u) -> Identity(),(t,u)-> I)
Out[39]:
The constraint operators for the first two equations are the usual ones for a stationary body and are precomputed. There are no constraints or constraint forces for the last three equations.
In [40]:
function TimeMarching.plan_constraints(u::TU,t::Float64)
B₁ᵀ, B₂ = TimeMarching.plan_constraints(u[1],t,sys) # These are used by both the first and second equations
return (B₁ᵀ,B₁ᵀ,f->zero(u[3]),f->zero(u[4]),f->zero(u[5])), (B₂,B₂,u->Vector{Float64}(),u->Vector{Float64}(),u->Vector{Float64}())
end
Set up the integrator here
In [41]:
@time ifherk = IFHERK(u,f,sys.Δt,plans,TimeMarching.plan_constraints,
(TimeMarching.r₁,TimeMarching.r₂),rk=TimeMarching.RK31,isstored=true)
Out[41]:
Initialize the state vector and the history vectors
In [57]:
t = 0.0
u = (zero(w₀),zero(w₀),zero(ΔX),zero(ΔX),[0.0,0.0])
thist = Float64[];
w1hist = [];
w2hist = [];
dX1hist = [];
dX2hist = [];
X̂hist = [];
tsample = Δt; # rate at which to store field data
Set the time range to integrate over.
In [58]:
tf = 8π
T = Δt:Δt:tf;
In [59]:
@time for ti in T
global t, u, f = ifherk(t,u)
push!(thist,t)
if (isapprox(mod(t,tsample),0,atol=1e-6) || isapprox(mod(t,tsample),tsample,atol=1e-6))
push!(w1hist,deepcopy(u[1]))
push!(w2hist,deepcopy(u[2]))
push!(dX1hist,deepcopy(u[3]))
push!(dX2hist,deepcopy(u[4]))
push!(X̂hist,deepcopy(u[5]))
end
end
println("solution completed through time t = ",t)
In [60]:
function vorticity(w::Nodes{Fields.Dual},sys)
ω = deepcopy(w)
return rmul!(ω,1/cellsize(sys))
end
function velocity(w::Nodes{Fields.Dual},sys)
return rmul!(curl(sys.L\w),-1)
end
function streamfunction(w::Nodes{Fields.Dual},sys)
return rmul!(sys.L\w,-cellsize(sys))
end
function nl(w::Nodes{Fields.Dual},sys)
return rmul!(TimeMarching.r₁(w,0.0,sys),1/cellsize(sys))
end
Out[60]:
In [61]:
function average(hist::Vector{Any},itr)
avg = zero(hist[1])
return avg .= mapreduce(x -> x/length(itr),+,hist[itr]);
end
Out[61]:
In [66]:
iplot = length(w1hist) # index of time step for plotting
ω₁ = vorticity(w1hist[iplot],sys)
q₁ = velocity(w1hist[iplot],sys)
ψ₁ = streamfunction(w1hist[iplot],sys)
nothing
In [67]:
xg,yg = coordinates(w₀,sys.grid)
plot(xg,yg,ω₁,levels=range(-10,10,length=30),color=:RdBu)
plot!(body)
Out[67]:
In [68]:
plot(xg,yg,ψ₁,levels=range(-1,1,length=31),color=:RdBu)
plot!(body,fill=:false)
Out[68]:
In [72]:
ix = 201
plot(yg,ω₁[ix,:],ylim=(-10,10),xlim=(1,4),label="DNS",xlabel="y")
plot!(yg,map(y -> ω₁ex(xg[ix],y,thist[iplot]),abs.(yg)),label="Analytical")
#plot!(yg,-real(W₁.(abs.(yg))*exp(-im*thist[iplot])))
Out[72]:
The following functions pick off the vorticity and velocity histories at an index pair. Note that u and v components that share the same indices are at different physical locations, and each are at different locations from vorticity.
In [73]:
function vorticity_history(sample_index::Tuple{Int,Int},whist,sys)
i, j = sample_index
return map(x -> vorticity(x,sys)[i,j],whist)
end
function velocity_history(sample_index::Tuple{Int,Int},whist,sys)
i, j = sample_index
tmp = reshape(collect(Iterators.flatten(
map(x -> (q = velocity(x,sys); [q.u[i,j], q.v[i,j]]),whist))),
2,length(whist))
return tmp[1,:], tmp[2,:]
end
# This computes the history of the rhs of the second-order equation at the sample point
function nl_history(sample_index::Tuple{Int,Int},whist,sys)
i, j = sample_index
return map(w1 -> TimeMarching.r₁(w1,0.0,sys)[i,j]/cellsize(sys),whist)
end
function centroid_history(Xhist)
tmp = reshape(collect(Iterators.flatten(Xhist)),2,length(Xhist))
return tmp[1,:], tmp[2,:]
end
Out[73]:
In [78]:
sampij = (201,135) # indices of sample point
Out[78]:
In [79]:
ω₁hist = vorticity_history(sampij,w1hist,sys);
Get coordinates of sample point on dual node grid
In [80]:
xg, yg = coordinates(w₀,sys.grid)
xeval = xg[sampij[1]]
yeval = yg[sampij[2]]
reval = sqrt(xeval^2+yeval^2)
coseval = xeval/reval
sineval = yeval/reval
reval
Out[80]:
In [81]:
plot(thist,ω₁hist,label="DNS",ylim=(-10,10))
plot!(thist,ω₁ex(xeval,yeval,thist),label="Analytical") # note that we correct the sign of W₁ here
#plot!(thist,real.(-W₁(reval)*sineval*exp.(-im*thist)),label="Analytical") # note that we correct the sign of W₁ here
Out[81]:
In [82]:
u₁hist, v₁hist = velocity_history(sampij,w1hist,sys);
Get the coordinates of the sample point in the u-component edges for plotting this component
In [83]:
xuedge, yuedge, xvedge, yvedge = coordinates(Edges(Primal,w₀),sys.grid)
xeval = xuedge[sampij[1]]
yeval = yuedge[sampij[2]]
reval = sqrt(xeval^2+yeval^2)
coseval = xeval/reval
sineval = yeval/reval
Out[83]:
In [84]:
plot(thist,u₁hist,ylim=(-2,2),label="DNS")
plot!(thist,u₁ex(xeval,yeval,thist),xlim=(0,8π),label="Analytical")
Out[84]:
Plot the centroid history
In [85]:
xhist, yhist = centroid_history(X̂hist)
plot(thist,xhist,label="X̂c")
plot!(thist,yhist,label="Ŷc")
Out[85]:
In [86]:
iplot = length(w2hist)
ω₂ = vorticity(w2hist[iplot],sys)
q₂ = velocity(w2hist[iplot],sys)
ψ₂ = streamfunction(w2hist[iplot],sys)
nothing
In [89]:
xg,yg = coordinates(w₀,sys.grid)
plot(xg,yg,ω₂,levels=range(-5,5,length=30),color=:RdBu,clim=(-5,5))
plot!(body)
Out[89]:
Compute the average vorticity field over a period, $\overline{w}_2$
In [90]:
iplot = length(w2hist)
itr = iplot-floor(Int,2π/(Ω*Δt))+1:iplot
w2avg = average(w2hist,itr)
ω̄₂ = vorticity(w2avg,sys)
ψ̄₂ = streamfunction(w2avg,sys);
In [91]:
ψ̄₂exfield = Nodes(Fields.Dual,w₀)
ψ̄₂exfield .= [ψ̄₂ex(x,y) for x in xg, y in yg];
In [94]:
plot(xg,yg,ψ̄₂,levels=range(-0.2,0.1,length=31),clim=(1,2),xlim=(0,4),ylim=(0,4))
#plot!(xg,yg,ψ̄₂exfield,levels=range(-0.2,0.1,length=31),clim=(-2,-2),xlim=(0,4),ylim=(0,4))
plot!(body,fillcolor=:black,linecolor=:black)
Out[94]:
In [95]:
sqrt(xg[119]^2+yg[119]^2)
Out[95]:
In [96]:
ig = 120
plot(sqrt.(xg[ig:end].^2+yg[ig:end].^2),map(i -> ψ̄₂[i,i],ig:length(xg)),ylim=(-1,1),xlim=(1,5),label="DNS",xlabel="r")
plot!(sqrt.(xg[ig:end].^2+yg[ig:end].^2),map((x,y) -> ψ̄₂ex(x,y),xg[ig:end],yg[ig:end]),label="Analytical",grid=:true)
Out[96]:
In [119]:
sampij = (240,240) # indices of sample point
Out[119]:
In [120]:
xg, yg = coordinates(w₀,sys.grid)
xeval = xg[sampij[1]]
yeval = yg[sampij[2]]
reval = sqrt(xeval^2+yeval^2)
Out[120]:
In [121]:
ω₂hist = vorticity_history(sampij,w2hist,sys);
In [124]:
plot(thist,ω₂hist,label="DNS",ylim=(-30,30))
plot!(thist,ω₂ex(xeval,yeval,thist),label="Analytical")
Out[124]:
In [125]:
xuedge, yuedge, xvedge, yvedge = coordinates(Edges(Primal,w₀),sys.grid)
xeval = xuedge[sampij[1]]
yeval = yuedge[sampij[2]]
reval = sqrt(xeval^2+yeval^2)
Out[125]:
In [126]:
u₂hist, v₂hist = velocity_history(sampij,w2hist,sys);
In [128]:
plot(thist,u₂hist,ylim=(-2,2),label="DNS")
plot!(thist,u₂ex(xeval,yeval,thist),xlim=(0,8π),label="Analytical")
Out[128]:
In [129]:
plot(thist,v₂hist,ylim=(-2,2),label="DNS")
plot!(thist,v₂ex(xeval,yeval,thist),xlim=(0,8π),label="Analytical")
Out[129]:
In [132]:
iplot = 380
println("Time/period = ",thist[iplot]/(2π))
usamp = (w1hist[iplot],w2hist[iplot],dX1hist[iplot],dX2hist[iplot],X̂hist[iplot])
rhs1samp = TimeMarching.r₁(usamp,thist[iplot])
rhs2samp = TimeMarching.r₂(usamp,thist[iplot]);
In [139]:
θ = range(0,2π,length=length(coords.u)+1)
plot(θ[1:end-1],rhs2samp[2].u,xlim=(0,2π),ylim=(-2,2),label="DNS",xlabel="θ")
plot!(θ[1:end-1],map((x,y) -> u₂ex(x,y,thist[iplot]),cos.(θ[1:end-1]),sin.(θ[1:end-1])),label="Analytical")
plot!(θ[1:end-1],rhs2samp[2].v,xlim=(0,2π),ylim=(-2,2),label="DNS",xlabel="θ")
plot!(θ[1:end-1],map((x,y) -> v₂ex(x,y,thist[iplot]),cos.(θ[1:end-1]),sin.(θ[1:end-1])),label="Analytical")
Out[139]:
In [142]:
sampij = (240,240) # indices of sample point
Out[142]:
In [143]:
xg, yg = coordinates(w₀,sys.grid)
xeval = xg[sampij[1]]
yeval = yg[sampij[2]]
reval = sqrt(xeval^2+yeval^2)
coseval = xeval/reval
sineval = yeval/reval
sin2eval = 2*sineval*coseval
reval
Out[143]:
In [144]:
rhshist = nl_history(sampij,w1hist,sys);
In [145]:
itr = length(thist)-floor(Int,2π/(Ω*Δt))+1:length(thist)
println("Mean value of DNS-computed RHS = ",Statistics.mean(rhshist[itr]))
println("Mean value analytical RHS = ",real(f₀(reval)/p.Re*sin2eval))
In [146]:
plot(thist,rhshist,label="DNS")
plot!(thist,real.(f₀(reval)/p.Re .+ g₀(reval)/p.Re*exp.(-2im*thist))*sin2eval,label="Analytical")
Out[146]:
In [147]:
itr = length(thist)-floor(Int,2π/(Ω*Δt))+1:length(thist)
vx = zero(w₀)
vy = zero(w₀)
Xx = zero(w₀)
Xy = zero(w₀)
ψd = zero(w₀)
ψdhist = []
for (i,dX) in enumerate(dX1hist)
cellshift!((vx,vy),curl(sys.L\w1hist[i])) # -v₁, on the dual nodes
cellshift!((Xx,Xy),dX) # ΔX₁, on the dual nodes
ψd .= rmul!(Xx∘vy-Xy∘vx,0.5) # 0.5(v₁ × ΔX₁)
push!(ψdhist,ψd)
end
ψd .= average(ψdhist,itr);
In [148]:
xg,yg = coordinates(w₀,sys.grid)
ψ̄₂d = deepcopy(ψ̄₂)
ψ̄₂d .+= ψd
plot(xg,yg,ψ̄₂d,levels=range(-0.2,0.1,length=31),clim=(1,2),xlim=(0,4),ylim=(0,4))
#plot!(xg,yg,ψ̄₂exfield,levels=range(-0.5,0.5,length=31),clim=(-2,-2),xlim=(0,4),ylim=(0,4))
plot!(body,fillcolor=:black,linecolor=:black)
Out[148]:
In [ ]: