force coefficients are half the correct values

[biao-geng-me]

First, thank you and congratulations on this amazing work! As someone who's been working with CFD on HPC using CPUs, the speed of WaterLily using GPUs is just jaw-dropping.

We are trying to incorprate WaterLily in our workflow. So far we've tested a circular cylinder case and a square cylinder case, both 3D. The circular case had accurate force coefficient compared to literature. However, when switched to the square shape, I had roughly half the correct force coefficients (reference). For example, the mean drag @ Re=100 is about 1.5, but I got ~0.75 from WaterLily. I calculated the coefficient using the built-in metrics function total_force

    force = WaterLily.total_force(sim)
    force./=(0.5sim.L*nz*sim.U^2) # scale the forces! nz is the dimension in the z axis

Confoundingly, this part was not changed when switching from circular to square cylinder.

I inspected the pressure field (shown below) at one time instant, and it looked correct and compared well with the pressure field using another CFD code.

The drag force coefficient history is

My code for this case is attached if it helps. Thank you for your time! square_cylinder_3D.zip

[weymouth]

Apologies, I was teaching in the last quarter, so I'm just catching up with some of these issues. Did you figure this out?

[biao-geng-me]

@weymouth @marinlauber Thank you for the responses. No, we didn't figure it out. We decided to proceed with the project using WaterLily nontheless, since the flow field looked correct and we don't need the forces at the moment. It'll be great if you could take the time to look into this case. I lowkey hope it's not some dumb/simple mistaken in my setup, but my apologies if that's the case.

@marinlauber nz indeed is the spanwise dimension of the cylinder. I've double checked that the shape is indeed a square cylinder (not a cube). The spanwise dimension passes through the whole domain in the z direction.

Thank you again for your work and your input!

[weymouth]

I certainly hope this is just a simple mistake in the setup ;-)

As @Marin mentioned, the pressure outside the box looks perfect, and you're off by a factor of two. It's probably a bug in the input or output scaling.

1. Are you sure that sim.L is the same width used in the literature to scale the forces? From you picture, the full width looks to be around 50 cells...
2. Are you sure you've matched your reference Re correctly? 50 cells is fairly low, so hopefully you're trying to match a low Re case.
3. One good way to catch scaling issues is to play around with the geometry and resolution. How does this change with refinement, or when you adjust the span to width ratiio?

[weymouth]

Here is a very minimal working example:

using WaterLily,StaticArrays,TypedTables,Plots
function make_square(N::NTuple{n,Int};H=N[2]÷8,Re=200) where n
  U = ntuple(i->i==1 ? 1 : 0,n)       # n-dimensional flow
  perdir = n==3 ? (3,) : () # optional
  function sdf(x,t)
    y = abs.(SA[x[1],x[2]] .- N[2]÷2) # positive quadrant vector from the center
    √sum(abs2,y-min.(y,H))             # using half-width to define shape from center
  end
  Simulation(N,U,2H;body=AutoBody(sdf), ν=2H/Re, perdir) # use full-width for scaling
end

sim = make_square((256,128)); # 2D
hist = map(100:0.1:110) do t
  sim_step!(sim,t,remeasure=false)
  f = WaterLily.total_force(sim)/(0.5sim.L*sim.U^2)
  (t=t,drag=-f[1],lift=f[2])
end |> Table
plot(hist.t,hist.drag,label="drag")
plot!(hist.t,hist.lift,label="lift")
plot!(ylabel="F/½ρU²H",xlabel="tU/H",
  title="2D square forces at Re=UH/ν=200")

Obviously, this is low resolution and low Re, but the forces seem to have the right magnitude. Is this helpful?

[weymouth]

Here's the 3D version at Re=2000