Unphysical forces for thin bodies

[dfloryan]

TL;DR
It seems that ParametricBodies struggles with thin bodies.

The issue

We've run into the following issue: for a flapping NACA foil, the lift force is unphysical. We get the same issue if the body is instead an ellipse, but only if the ellipse is thin enough.

Representative result

Here's the lift on a NACA0010 foil. The foil is parallel to the freestream, sinusoidally heaving with a reduced frequency of 5 and amplitude equal to 5% of the foil's length. The Reynolds number is 100, and the length of the foil is 64 grid points.

Observation: the "corners" in the lift curve coincide with the trailing edge moving into a different cell (i.e., they coincide with changes in floor(heave)). Other corners coincide with changes in round(heave).

Here's the same result for a higher resolution (length of foil = 128 grid points). The corners in the lift curve again coincide with the trailing edge moving into a different cell, which occurs more frequently due to the more refined grid.

Notes

• This issue is independent of the Reynolds number.
• For Re = 100, a grid resolution of 64 points per chord should be sufficient to resolve the flow. Even though refining the grid makes the unphysical spikes smaller, I don't think the solution is to simply continue refining the grid until the spikes are gone. Even at a resolution of 192 grid points per chord—which is far above what is typically required for such a problem at a Reynolds number of 100—the lift curve still has unphysical spikes.
• Looking into the force calculation, the output of WaterLily.nds looks reasonable. The issue seems to be with the pressure field.
• The velocity and vorticity fields look reasonable despite the unphysical lift.
• The same issue occurs for a very thin ellipse (e.g., aspect ratio of 10). Thicker ellipses (e.g., aspect ratio of 5) don't have this issue.

Thoughts and minimal working example

We aren't sure if this is a WaterLily issue or a ParametricBodies issue. However, the following test suggests it may be a ParametricBodies issue.

Here's an example identical to the NACA example up top, except the body is a plate with some thickness (same as in the TwoD_Hover example). This allows us to make the body in two ways: using an SDF, and using ParametricBodies. We've found that for most values of the thickness thk, the two ways of creating the body give nearly identical lift curves. However, for the default value in the script below, the body created using ParametricBodies produces unphysical lift, while the one created with an SDF is physical. The two bodies look the same when we plot them.

using WaterLily,StaticArrays
using Plots
using ParametricBodies

# Calculate force due to pressure
function pressure_force(sim;t=WaterLily.time(sim.flow),T=promote_type(Float64,eltype(sim.flow.p)))
   sim.flow.f .= zero(eltype(sim.flow.p))
   @WaterLily.loop sim.flow.f[I,:] .= sim.flow.p[I]*WaterLily.nds(sim.body,loc(0,I,T),t) over I ∈ WaterLily.inside(sim.flow.p)
   sum(T,sim.flow.f,dims=ntuple(i->i,ndims(sim.flow.p)))[:] |> Array
end

# Simulation using SDF
function make_sim(;L=64,Re=1e2,sigma_L=1,U=1,thk=5+√2,n=3,m=2,T=Float64,mem=Array)
    nose = SA[L,L]
    h₀=T(0.05*L); 
    ω = T(2*sigma_L/L)
    # Line segment SDF
    function sdf(x,t)
        y = x .- SA[clamp(x[1],0,L),0]
        √sum(abs2,y)-thk/2
    end
    # Oscillating motion
    function map(x,t)
        h = SA[0, h₀ * sin(ω * t)]
        (x - nose - h)
    end
    Simulation((n*L,m*L),(U,0),L;ν=U*L/Re,body=AutoBody(sdf,map),T,mem)
end

# Simulation using ParametricBodies
function make_sim1(;L=64,Re=1e2,sigma_L=1,U=1,thk=5+√2,n=3,m=2,T=Float64,mem=Array)

    # Map from simulation coordinate x to surface coordinate ξ
    nose = SA[3*L/2,L]
    h₀=T(0.05*L); 
    ω = T(2*sigma_L/L)
    function map(x, t)
        h = SA[0, h₀ * sin(ω * t)]
        ξ = x - nose - h
        return SA[ξ[1], abs(ξ[2])]
    end

    # Define foil using NACA0010 profile equation: https://tinyurl.com/NACA00xx
    foil(s,t) = (s<π/4)*SA[L/2 + thk/2*cos(2*s),thk/2*sin(2*s)] + (s>3*π/4)*SA[-L/2 + thk/2*cos(2*s-π),thk/2*sin(2*s-π)] + (s>=π/4 && s<=3*π/4)*SA[L/2 - L*(s - π/4)*2/π,thk/2]
    body = HashedBody(foil,(0,π);map,T,mem)
    
    Simulation((n*L,m*L),(U,0),L;ν=U*L/Re,body,T,mem);
end

# Run and compare the two simulations
res=64;
sigma=5;
thk=5+√2
sim = make_sim(L=res,sigma_L=sigma,thk=thk);
sim1 = make_sim1(L=res,sigma_L=sigma,thk=thk);
nosc=1; tmax=(π/(sigma))*nosc

cycle = range(0,tmax,nosc*50)

p,p1 = [],[]

for t ∈ sim_time(sim).+cycle
    sim_step!(sim,t,remeasure=true)
    push!(p,pressure_force(sim))
end

for t ∈ sim_time(sim1).+cycle
    sim_step!(sim1,t,remeasure=true)
    push!(p1,pressure_force(sim1))
end

plot(cycle, 2last.(p)/sim.L, xlabel="tU/L", ylabel="Pressure lift")
plot!(cycle, 2last.(p1)/sim1.L, xlabel="tU/L", ylabel="Pressure lift")

[weymouth]

Thanks for the careful write-up.

First, the force routines in ParametricBodies were causing us some issues as well, so we've deprecated those functions. If you are using the current master branch of ParametricBodies, the force routines should be those from WaterLily.

In the bottom example it looks like you have removed that potential difference by writing something similar to the WaterLily function. So the only source of difference should be the body. I think I'll play around with this example and try to reduce it even further.

[weymouth]

I did this as an initial check. As expected, using the perfect locator gives identical results regardless of the body:

using WaterLily,StaticArrays
using Plots
using ParametricBodies

# Simulation
function make_sim(parametric=true;L=64,U=1,thk=3+√2,T=Float64,mem=Array)
    # Oscillating motion
    nose,h₀ = SA[L,L],T(0.05L)
    map(x,t) = x-nose-SA[0, h₀*sin(t*U/L)]

    # Line segment from [0,L] with thickness thk
    function sdf(x,t)
        y = x-SA[clamp(x[1],0,L),0]
        √(y'y)-thk/2
    end
    line(u,t) = SA[u,0]
    locate(x::SVector{2},t) = clamp(x[1],0,L)
    # gross hack you don't need when using Hashed or NURBS bodies
    @eval ParametricBodies.notC¹(::Function,u) = true 

    # Define bodies
    auto_body = AutoBody(sdf,map)
    para_body = ParametricBody(line,locate;map,thk,boundary=false) 

    Simulation((3L,2L), (U,0), L; ν=U*L/1e2, T, mem, body = parametric ? para_body : auto_body)
end
function get_lift!(sim,t)
    sim_step!(sim,t)
    force = WaterLily.pressure_force(sim) # make sure you are using ParametricBodies#master
    force[2]/(0.5sim.L*sim.U^2) # lift coefficient
end

# Run and compare the two simulations
time = 0.1:0.1:20 # time scale is sim.L/sim.U
plt = plot(xlabel="tU/L",ylabel="Pressure lift");
for parametric in (true,false)
    # Create simulation 
    sim = make_sim(parametric,L=32)

    # Simulate through the time range and get forces
    lift = [get_lift!(sim,t) for t in time]

    # Plot
    plot!(time, lift, label="Parametric=$parametric")
end; plt

Now I'll check out the other body & locator types.

[weymouth]

Replacing the parametric body definition with

line(u,t) = SA[u,0]
para_body = HashedBody(line,(0,L);map,thk,boundary=false) 

gives an identical result.

[weymouth]

line = BSplineCurve(SA{T}[0 L; 0 0]; degree=1)
para_body = ParametricBody(line;map,thk,boundary=false)

is also perfect. Could you try to start from one of these examples and morph it towards what you need to isolate what's breaking?

[dfloryan]

Thanks, I get the same results as you. I modified your script to use a NACA0010:

using WaterLily,StaticArrays
using Plots
using ParametricBodies

# Simulation
function make_sim(;L=64,U=1,T=Float64,mem=Array)
    # Oscillating motion
    nose,h₀ = SA[L,L],T(0.05L)
    map(x,t) = x-nose-SA[0, h₀*sin(t*U/L)]
    
    # NACA0010 foil
    NACA(s) = 0.5f0*(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)]
    para_body = HashedBody(foil,(0,1);map,boundary=false)

    Simulation((3L,2L), (U,0), L; ν=U*L/1e2, T, mem, body = para_body)
end
function get_lift!(sim,t)
    sim_step!(sim,t)
    force = WaterLily.pressure_force(sim) # make sure you are using ParametricBodies#master
    force[2]/(0.5sim.L*sim.U^2) # lift coefficient
end

# Run simulation
time = 0.1:0.1:10 # time scale is sim.L/sim.U
plt = plot(xlabel="tU/L",ylabel="Pressure lift");

sim = make_sim(L=32)

# Simulate through the time range and get forces
lift = [get_lift!(sim,t) for t in time]

# Plot
plot!(time, lift)
plt

which produces the following lift:

Definitely unphysical. I also tried removing the boundary=False argument from HashedBody since I hadn't seen it before. That made the simulation an order of magnitude slower and produced unphysically large lift (see below), which suggests the pressure solver had issues converging.

[weymouth]

Since the issue is with the specific body shape, it will be helpful to look directly at the SDF and not just at the forces. When we look at the SDF of my line body using this code

sim = make_sim(L=32)
measure_sdf!(sim.flow.σ,sim.body)
contour(sim.flow.σ',levels=[-sim.ϵ,0,sim.ϵ], clims = (-2sim.ϵ,2sim.ϵ))

we see

where, critically, there is a sdf≤-ϵ region within the body. This region is required for BDIM-σ to enforce the pressure boundary condition. In this image I've set thk=2ϵ+√2 (which is the minimum for the worst case placement and orientation), and use the WaterLily default ϵ=1 for the BDIM convolution thickness. See Marin's paper https://arxiv.org/abs/2110.06535 where he does a great job developing BDIM-σ for thin bodies.

So the question is how to achieve a thin sectional shape? First, Marin gave himself a little more room by setting ϵ=1/2 in his paper. This only buys you one grid cell with BDIM-σ, but every little bit helps. Then you can use a variable thickness on either side of a defined centerline. You need to add ParametricBodies#25 for that since PR #25 hasn't been merged yet (but will be soon, I hope). Here's the code:

function make_sim(;L=64,U=1,T=Float64,ϵ=T(0.5),thk=2ϵ+√2,mem=Array)
    # Oscillating motion
    nose,h₀ = SA[L,L],T(0.05L)
    map(x,t) = x-nose-SA[0, h₀*sin(t*U/L)]

    # NACA0010 foil
    NACA(s) = L*(0.2969f0s-0.126f0s^2-0.3516f0s^4+0.2843f0s^6-0.1036f0s^8) # thickness (not half)
    centerline(s,t) = SA[L*(1-s)^2,0]
    para_body = para_body = HashedBody(centerline,(0,1);map,thk=NACA,boundary=false)

    Simulation((3L,2L), (U,0), L; ν=U*L/1e2, T, mem, body = para_body, ϵ)
end

and here's the SDF:

Notice that a bunch of the trailing edge is cut off since the thickness is less than the minimum. Marin suggested to use max(thk,NACA), but that basically turns the back of the body into the line segment.

function NACA(s) 
        t = L*(0.2969f0s-0.126f0s^2-0.3516f0s^4+0.2843f0s^6-0.1036f0s^8) # thickness (not half)
        s>0.5 ? t : max(thk,t) # ensure min thickness at tail
end

And gives this as the lift:

Probably not ideal, but safe.

Audrey's paper (https://eprints.soton.ac.uk/369635/1/Maertens2014%2520CMAME.pdf) dealt with the thin trailing edge section a different way. By subtracting the SDFs defined by the top and bottom contours of the section, and dealing with the BDIM kernels carefully (see Appendix A), she managed to get really nice trailing edge behaviour. This hasn't been implemented in WaterLily yet, but it could be an interesting extension.

[dfloryan]

Thanks, that makes sense.

For completeness, I pulled the latest commit, which includes the PR you mentioned. I was able to reproduce the second contour plot you made (the SDF of the NACA0010 foil with its trailing edge cut off). Below is the lift for the foil with the cut off trailing edge:

And here it is when the resolution is doubled (L = 64):

The lift has non-physical jumps that become smaller and more frequent as the resolution is increased. I'm not sure whether this behaviour is to be expected, given how the SDF looks by the trailing edge.

In any case, this SDF does a poor job of representing the true foil because of how thin it is near the trailing edge, so it may not be the right approach anyway. It seems that the approach described in Audrey's paper is the right one for bodies with sharp corners. It would be awesome to see it implemented in WaterLily!

FYI, here's the script I used to produce the plots:

using WaterLily,StaticArrays
using Plots
using ParametricBodies

# Simulation
function make_sim(;L=64,U=1,T=Float64,ϵ=T(0.5),thk=2ϵ+√2,mem=Array)
    # Oscillating motion
    nose,h₀ = SA[L,L],T(0.05L)
    map(x,t) = x-nose-SA[0, h₀*sin(t*U/L)]

    # NACA0010 foil
    NACA(s) = L*(0.2969f0s-0.126f0s^2-0.3516f0s^4+0.2843f0s^6-0.1036f0s^8) # thickness (not half)
    centerline(s,t) = SA[L*(1-s)^2,0]
    para_body = HashedBody(centerline,(0,1);map,thk=NACA,boundary=false)

    Simulation((3L,2L), (U,0), L; ν=U*L/1e2, T, mem, body = para_body, ϵ)
end
function get_lift!(sim,t)
    sim_step!(sim,t)
    force = WaterLily.pressure_force(sim) # make sure you are using ParametricBodies#master
    force[2]/(0.5sim.L*sim.U^2) # lift coefficient
end

# Run simulation
time = 0.01:0.01:10 # time scale is sim.L/sim.U
plt = plot(xlabel="tU/L",ylabel="Pressure lift");

sim = make_sim(L=32)
# Simulate through the time range and get forces
lift = [get_lift!(sim,t) for t in time]

# Plot
plot!(time, lift)
plt

[weymouth]

Yes, I think your choices are to either use the minimal thickening that Marin developed in his paper or the modified kernel moment functions the Audrey developed in hers. Ignoring the issue doesn't work, as you've shown.

Audrey's version requires defining both surfaces separately in a new DifferenceBody struct, and then writing a new measure!(a::Flow,body::DifferencetBody) function to use the modified moments. The data for a DifferenceBody could just be a tuple with two AbstractBodies. This is pretty close to the Bodies types. https://github.com/WaterLily-jl/WaterLily.jl/blob/080d08d7e8a01b80ea21248f1c31090cacc8b4fa/src/AutoBody.jl#L55

Why don't you open an issue on WaterLily in this direction and reference this issue. The suggested change will be there, not here.