WaterLily-jl / WaterLily-Examples

WaterLily tutorial files
GNU General Public License v3.0
13 stars 5 forks source link

How to get the boundary pressure value on the left foil? #3

Open liangaomng opened 2 weeks ago

liangaomng commented 2 weeks ago
input parameter in line 8

using WaterLily,StaticArrays using ParametricBodies # New package using DelimitedFiles function make_sim(;L=48,Re=1e3,St=0.3,αₘ=-π/18,U=1,T=Float32,mem=Array)

Map from simulation coordinate x to surface coordinate ξ

nose,pivot = SA[L,2L],SA[0.25f0L,0]

######input parameter
θ₀ = T(8/180*π);h₀=T(L);ω=0*T(π*St*U/h₀);vy_front=T(U/L/10);dis=T(3L)#duration=1
######
function map(x,t)
    back = x[1]>nose[1]+L       # back body?

    h = back ? 0 : h₀*t*vy_front  # y shift
    S = back ? 3L : 3L-dis           # horizontal shift
    θ = back ? θ₀*cos(ω*t) : T(30/180*π) 
    R = SA[cos(θ) -sin(θ); sin(θ) cos(θ)]
    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)
#test

Simulation((8L, 4L), (U, 0), L; ν = U * L / Re, body, T, mem)

#Simulation((8L,4L),(U,0),L;ν=U*L/Re,body,T,mem)

end

using CUDA

using Plots; gr()

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

addbody(x,y;c=:black) = Plots.plot!(Shape(x,y), c=c, legend=false) 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) end

function sim_gif!(sim;duration=1,step=0.1,verbose=true,R=inside(sim.flow.p), remeasure=false,plotbody=false,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...) plotbody && body_plot!(sim)

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

    plot_boundary(sim.body,tᵢ)

     # Save flow.σ and flow.p to text files--lam
    # open("flow_sigma_$(round(tᵢ, digits=2)).txt", "w") do file
    # end
    # open("flow_p_$(round(tᵢ, digits=2)).txt", "w") do file
    #     writedlm(file, sim.flow.p |> Array)

    # end

end

end

sim_gif!(sim,duration=0.2,step=0.1,clims=(-16,16),remeasure=true, plotbody=true,shift=(-2,-1.5),axis=([], false),cfill=:seismic, legend=false,border=:none)

Hi,Porf , very nice waterlily. attached is my code, I want to get the motion like coord on the left foil[x,y,t] and right foil' coord.

b-fg commented 1 week ago

To use multiple ParametricBodis and be able to compute their unique forces you should use https://github.com/WaterLily-jl/WaterLily.jl/pull/158. This is still a draft, so not properly tested for production.

liangaomng commented 1 week ago

Hi Prof, Thanks for your reply. I have another question: if we consider a single foil, how can we extract the pressure values along its boundary? These should be scalar values accompanied by their respective x and y coordinates.

b-fg commented 1 week ago

Take a look at the pressure force routine:

function pressure_force(p,df,body,t=0,T=promote_type(Float64,eltype(p)))
    df .= zero(eltype(p))
    @loop df[I,:] .= p[I]*nds(body,loc(0,I,T),t) over I ∈ inside(p)
    sum(T,df,dims=ntuple(i->i,ndims(p)))[:] |> Array
end

You just need to not integrate the force, hence look at df values, and get the coordinates for when df!=0 using loc.

liangaomng commented 1 week ago

Thanks for your help! If for a multibody class, how can I separate and retrieve the forces for each body? Do you have an example case that I could directly apply to my foil code?

b-fg commented 1 week ago

For regular Bodies (non parametric), take a look at this example. Otherwise, for multiple parametric bodies, you can test that PR I linked before.