Validation of the flow past cylinder case (re=100,200)

[liangaomng]

Hi Prof, Thanks for your effort. But I am confused that in the flow past cylinder case, I test the different grids, the waterlily calculate the mean of Cd is almost 1.53, and the rms cl is almost 0.32. But the numercial mean of Cd should be 1.35 and rms of cl should be 0.22. I do not know if there is problem? The validtion result is from "Quantitative numerical analysis of flow past a circular cylinder at Reynolds number between 50 and 200"

[b-fg]

Please share a reproducible example so we can look into it. Thanks!

[liangaomng]

Hi prof,

Attached is my code.

using WaterLily
using Plots
function circle(n,m;Re=2000,U=1)
    # signed distance function to circle
    radius, center = m/10, m/2-1
    sdf(x,t) = √sum(abs2, x .- center) - radius

    Simulation((n,m),   # domain size
               (U,0),   # domain velocity (& velocity scale)
               2radius; # length scale
               ν=U*2radius/Re,     # fluid viscosity
               body=AutoBody(sdf)) # geometry
end
function get_forces!(sim,t)
   sim_step!(sim,t,remeasure=false)
   println(t)
   force = WaterLily.pressure_force(sim)+WaterLily.viscous_force(sim)
   force./(0.5sim.L*sim.U^2) # scale the forces!
end
circ = circle(1024,512)
sim_step!(circ)
# Simulate through the time range and get forces
time = 1:0.05:80 # time scale is sim.L/sim.U
forces = [get_forces!(circ,t) for t in time];

#Plot it
plot(time,[first.(forces) last.(forces)],
   labels=["drag" "lift"],
   xlabel="tU/L",
   ylabel="Pressure force coefficients")

[b-fg]

(I just formatted the code for you, but please use triple quote blocks and Julia tag ```julia to help on readability).

[b-fg]

Just took a look. Here are a few comments:

• You mention that you use this to validate for "Reynolds number between 50 and 200", but the example you shared is using Re=2000.
• You should provide a link to the reference for the quantitative result you mention.
• The domain seems a bit small, I would inflate it +50% to be on the safe side.
• The example does not contain how you take the mean of the forces. Remember that when taking an average of the force (or any signal in general), you must discard the transient part.
• Also from the forces plot, it seems that only 13 wake cycles are captured. You need more than that for an accurate average.

Closing this now. Feel free to reopen if after all this the issue persists.

[b-fg]

For reference, here is the fixed version of this test case, with a bigger domain and at the correct Reynolds (Re=100). We get very close values to the reference data (Mean CD = -1.36; RMS CL = 0.24).

using Revise, WaterLily, Plots, CUDA, Distributions, Printf

function circle(n,m;Re=150,U=1,mem=Array,T=Float32)
    # signed distance function to circle
    radius, center = m/16, m/2
    sdf(x,t) = √sum(abs2, x .- center) - radius

    Simulation((n,m),   # domain size
               (U,0),   # domain velocity (& velocity scale)
               2radius; # length scale
               ν=U*2radius/Re,     # fluid viscosity
               body=AutoBody(sdf),  # geometry
               mem,
               T)
end

function get_forces!(sim,t)
   sim_step!(sim,t,remeasure=false)
   println(t)
   force = WaterLily.total_force(sim)
   force./(0.5sim.L*sim.U^2) # scale the forces!
end

n, m = 1024, 512
T = Float32
mem = CuArray
Re = 100
U = 1
sim = circle(n,m; U, Re, T, mem)
sim_step!(sim)

noise = rand(Uniform(-0.1*U,0.1*U), size(sim.flow.u)...) |> mem
sim.flow.u .+= T.(noise) # poke flow to speed up wake symmetry break

# Simulate through the time range and get forces
time = 1:0.05:125 # time scale is sim.L/sim.U
forces = [get_forces!(sim,t) for t in time]

# Plot it
Plots.plot(time,[first.(forces) last.(forces)],
   labels=["drag" "lift"],
   xlabel="tU/L",
   ylabel="Total force coefficients")

println("Mean CD = $(@sprintf("%2.2f", mean(first.(forces)[time.>75])))")
println("RMS CL = $(@sprintf("%2.2f", std(last.(forces)[time.>75])))")

Note that, at this Re, you can also run it with half the resolution and get still good results (Mean CD = -1.34; RMS CL = 0.24)

[liangaomng]

hi professor,

Thanks so much. BUT I still have problem，when I do the simulation on naca0012， I can simulate the vortex shedding.

When I extract force, every time I extract it, it is different, and even nan appears. Do you know why?

using WaterLily
using StaticArrays
using Plots
using WaterLily,StaticArrays
using ParametricBodies
using DelimitedFiles
using Plots

using Plots; gr()
function body_plot!(sim;levels=[0],lines=:black,R=inside(sim.flow.p))
   WaterLily.measure_sdf!(sim.flow.σ,sim.body,WaterLily.time(sim))
   contour!(sim.flow.σ[R]'|>Array;levels,lines)

   # body 
   data = sim.flow.σ[R]' |> Array

end

function flood(f::Array;shift=(0.,0.),cfill=:RdBu_11,clims=(),levels=10,kv...)
    if length(clims)==2
        @assert clims[1]<clims[2]
        @. f=min(clims[2],max(clims[1],f))
    else
        clims = (minimum(f),maximum(f))
    end
    Plots.contourf(axes(f,1).+shift[1],axes(f,2).+shift[2],f',
        linewidth=0, levels=levels, color=cfill, clims = clims, 
        aspect_ratio=:equal; kv...)
end

function make_sim(;L=48,Re=1e3,St=0.3,U=1,T=Float32,mem=Array,angles=[0,0])
   # Map from simulation coordinate x to surface coordinate ξ
   nose,pivot = SA[1L,3.5L],SA[0.f0L,0]
   #前后攻角
   fist_aoa = angles[1]
  
   ######input parameter
   h₀=T(0.5L);ω=T(π*St*U/h₀);vy_front=0;#duration=1
   Dis = 2L #相距
   Ini = 0L
   ###得到全局坐标
   function map(x,t)

       back = false      # back body?
       h = back ? 0 : Ini+h₀*t*vy_front  # y shift
       S = back ? Dis : 0           # horizontal shift
       θ1 = back ? T(0/180*π) : T(fist_aoa/180*π) 
       R = SA[cos(θ1) -sin(θ1); sin(θ1) cos(θ1)]
       h = SA[S,h]#h₀*sin(ω*t)]
       ξ = R*(x-nose-h-pivot)+pivot # move to origin and align with x-axis
       return SA[ξ[1],abs(ξ[2])]    # reflect to positive y
   end

   # Define foil using NACA0012 profile equation: https://tinyurl.com/NACA00xx
   NACA(s) = 0.6f0*(0.2969f0s-0.126f0s^2-0.3516f0s^4+0.2843f0s^6-0.1036f0s^8)
   foil(s,t) = L*SA[(1-s)^2,NACA(1-s)]
   

   body = HashedBody(foil,(0,1);map,T,mem)


   Simulation((8L,7L),(U,0),L;ν=U*L/Re,body,T,mem)
end
function sim_gif!(sim;duration=1,step=0.1,verbose=true,R=inside(sim.flow.p),
   remeasure=true,plotbody=true,kv...)
   t₀ = round(sim_time(sim))
   @time @gif for tᵢ in range(t₀,t₀+duration;step)
      sim_step!(sim,tᵢ;remeasure)
      @inside sim.flow.σ[I] = WaterLily.curl(3,I,sim.flow.u)*sim.L/sim.U
      flood(sim.flow.σ[R]|>Array; kv...)
      body_plot!(sim)

      verbose && println("tU/L=",round(tᵢ,digits=4),
      ", Δt=",round(sim.flow.Δt[end],digits=3))

   end


end
# intialize and run
sim = make_sim()
#sim_gif!(sim,duration=1,clims=(-10,10),plotbody=true)
function get_forces!(sim,t)
   println("t",t)
   sim_step!(sim,t,remeasure=false)
   force = WaterLily.viscous_force(sim.flow,sim.body)+WaterLily.pressure_force(sim.flow,sim.body)
   println("t",force)
   force./(0.5sim.L*sim.U^2) # scale the forces!
end

# Simulate through the time range and get forces
time = 1:0.1:50 # time scale is sim.L/sim.U
forces = [get_forces!(sim,t) for t in time];

#Plot it
plot(time,[first.(forces) last.(forces)],
   labels=["drag" "lift"],
   xlabel="tU/L",
   ylabel="Pressure force coefficients")

[b-fg]

This is a different issue. So please create a new one if needed.

[liangaomng]

Get it! Thanks