Erratic Pressure Behaviour

[brmcatee]

I'm experiencing erratic (not physically possible) pressure fluctuations throughout the simulation of flapping foils (almost identical to the example provided on GitHub). The issue appears to begin happening at a certain grid size (i.e., works well at small grid sizes of say L = 32, (x,y) = (8L,4L), but no longer works at, for example, L = 64, (x,y) = (16L,10L)). Velocity values are unaffected, indicating an issue with the pressure solver. Here is the code for my minimal working example:

using WaterLily,StaticArrays, ParametricBodies, CSV, DataFrames, Plots; gr()

function make_sim(;L=64,Re=7e3,St=0.35,αₘ=-π/18,U=1,T=Float32,mem=Array)
    # Map from simulation coordinate x to surface coordinate ξ
    nose,pivot = SA[6L,5L],SA[0.25f0L,0] # nose starting coordinates, pivot point (1/4 chord)
    θ₀ = T(αₘ+atan(π*St)); h₀=T(L); ω=T(π*St*U/h₀) 
    function map(x,t)
        θ = θ₀*cos(ω*t); R = [cos(θ) -sin(θ); sin(θ) cos(θ)]
        h = SA[0,h₀*sin(ω*t)]
        ξ = R*(x-nose-h-pivot)+pivot
        return SA[ξ[1],abs(ξ[2])]
    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((16L,10L),(U,0),L;ν=U*L/Re,body,exitBC=true,T,mem) # Note exitBC = true
end

# Iniialise the simulation
sim = make_sim()

# Create an array to store pressure values
p_field = []

# Step through the simulation
time = 0:0.1:5pi
for t in time
    sim_step!(sim, t, remeasure = true)
    pressure_field = sim.flow.p
    # Choosing a location near the inlet as 'free stream' pressure
    push!(p_field, pressure_field[20,20])
end

# Export data in a CSV file
df_pinf = DataFrame(value = p_field)
CSV.write("pressure_test.csv", df_pinf)

The data from the .csv file produced from the above code results in the following plot: This is the pressure at a fixed location near the inlet.

Other observations:

• The issue occurs regardless of whether the domain size is a power of 2 or not
• The magnitude of the fluctuations generally increase with domain size

[gusjdn]

Please see issue 36 : https://github.com/WaterLily-jl/WaterLily.jl/issues/36 We experienced a similar problem. Maybe Gabe can give more details, he made some changes in the BCs which helped, but the problem did not disappear. The issue was not defined as a priority back then. In fact, we experienced a similar problem just a few days ago.

Have you tried changing the exitBC setting?

[brmcatee]

I've experienced the issue both when exitBC is true and false.

[brmcatee]

I'm experiencing the same drift as you did in #36 , it just isn't visible in the above figure due to the scale of the drift versus the scale of the major fluctuations. The following figure shows a similar phenomenon as you outlined. The different curves represent different (x,y) sizes, ranging from (16L,8L) to (64L,16L) (and also the placements of the foil in the domain), where the rate of pressure loss generally decreases with an increase in domain size, and this trend makes sense.

Furthermore, the rate of pressure decrease does not seem to be linear: The above figure shows two different domain sizes (the data has been averaged using movmean in matlab to better observe the trend). Increasing the domain size from any of the above curves typically triggers the major fluctuations to begin occurring.

[weymouth]

Thanks for bringing this up. In that link referenced above #36, I explained that the value of pressure isn't uniquely defined up to a scalar offset in these simulation. The absolute value of pressure has no influence on the velocity or forces.

That doesn't mean everything is ok, it just means that you can't take a pressure reading, see wiggles, and know there is a problem.

If you pick one value of pressure, say p[2,2] and use it to define a pressure coefficient c_p = (p-p[2,2])/(0.5U^2), is there still a problem with drift or unphysical wiggles? If so, can you write a minimum reproducible example and we can try to figure out how to fix it?

[weymouth]

I should also mention that I'm in the process of merging in an update to the pressure solver @marinlauber worked on #135. He was trying to deal with some instabilities in the forces, and also give us better tools for diagnosing pressure problems. You could check that out.

[samuelsuntree]

Thanks for bringing this up. In that link referenced above #36, I explained that the value of pressure isn't uniquely defined up to a scalar offset in these simulation. The absolute value of pressure has no influence on the velocity or forces.

That doesn't mean everything is ok, it just means that you can't take a pressure reading, see wiggles, and know there is a problem.

If you pick one value of pressure, say p[2,2] and use it to define a pressure coefficient c_p = (p-p[2,2])/(0.5U^2), is there still a problem with drift or unphysical wiggles? If so, can you write a minimum reproducible example and we can try to figure out how to fix it?

Hi Gabe I am Boai :). I have been working with Marin on the pressure field for the past weeks. Actually calculating the integral of a convex body will remove the global shift so its good, but part of the output is to monitor a single pressure point in the field(like p[2,2]). Thus, for a quick fix on my side, is it a temporary solution that we calculate the mean pressure of the 4 corners(considering farest in the field) and remove the mean for every time step?

[weymouth]

Sure. just define p_inf = 0.5p[2,2]+0.5p[2,end-1] (or whatever) and then define c_p = (p .- p_inf)/(0.5sim.U^2). That should look pretty stable.