IB volume-integral force is wrong in the streamwise component for a thin inclined plate

[sbryngelson]

Summary

The immersed-boundary force that MFC reports (D/ib<id>_forces.dat, and the force/torque in the IB state files)
is computed as a volume integral of -grad p + div tau over the cells marked as body, in s_compute_ib_forces.
For a thin plate inclined to the mesh, the streamwise component of that integral is wrong by an order of
magnitude, and can come out with the wrong sign. The transverse component is fine.

Evidence

Case: rectangular flat plate, aspect ratio 4, pitching about its leading edge from 0 to 45 degrees on the
AIAA canonical smoothed ramp (K = pi/8), Re_c = 300, Ma 0.2, half span with a symmetry plane at the root,
WENO5 + HLLC, fd_order = 4.

A plate held at 45 degrees carries a force essentially normal to its surface, so C_D / C_L should be about 1.
It is, at the instant the rotation stops, and then the streamwise component collapses while the transverse one
holds:

time after ramp start	C_D/C_L, 2 cells across the section	4 cells across
1.0 (ramp end)	1.10	1.04
2.0	0.15	0.11
4.0	-0.22	0.10

A negative drag on a plate at 45 degrees is not physically possible.

Cross-checked with a control-volume momentum balance over a box enclosing the body,

F = closed integral over S of ( -p n + tau.n - rho u (u.n) ) dS  -  d/dt of the momentum inside S

evaluated entirely on grid planes away from the body, so it does not depend on the IB reconstruction. On a
static cylinder at Re 40 this balance returns C_D = 1.534 to 1.540 against the published 1.540 (Canuto & Taira
2015), insensitive to box size, so it is trustworthy. On the pitching plate, at times where the unsteady term is
small:

resolution	t	C_D from ib_forces	C_D from the control volume	C_L from ib_forces	C_L from the control volume
4 cells across	5.00	0.12	1.04	1.02	0.95
2 cells across	7.08	-0.22	0.85	0.86	0.87
2 cells across	8.33	-0.22	1.13	0.92	0.71

So the flow is approximately right and the reported streamwise force is not. Refining from 2 to 4 cells across
the section removes the negative sign but leaves the magnitude an order of magnitude low.

What I think is going on

I have not pinned the mechanism, so treat this as a report rather than a diagnosis. The volume integral equals
the surface traction only if the field inside the body is a smooth extension of the wall state. For a body a few
cells thick, every interior cell is a ghost point whose image point is reflected across the nearest surface, and
the fd_order = 4 stencil used in the integral reaches two cells beyond the body on both sides. My guess is
that the telescoping that should leave a surface term is incomplete in the direction across which the body is
thin, and at incidence that direction has a streamwise projection.

Things ruled out along the way: the wing root sitting on the symmetry plane (extending the plate through the
plane changed nothing), domain blockage (a wider domain moved results the other way), and grid resolution alone
(the error survives a 2x refinement).

Why it matters

Thin lifting surfaces are a natural use of the IB, and the force is usually the quantity of interest. Anyone
computing drag or thrust on a plate, wing or flapping surface with the current diagnostic will get a number that
looks plausible and is wrong.

Reproducer

examples/3D_ibm_pitchup_plate (added in #1847) at its default settings; the C_D/C_L ratio after the ramp is
the tell, and needs no reference data to see.

[sbryngelson]

Mechanism, now pinned

I said above that I had not pinned the mechanism. I have now, by reproducing the force from a saved field and
taking it apart. The flow is right; the ghost-filled interior pressure is wrong, and the volume integral is
built entirely out of it.

Same case as above, AR 4 plate held at 45 degrees, four cells across the section, Re 300, at t - t0 = 2.
All three numbers below come from the same post-processed field:

what is integrated	Fx	Fz	|F|
surface integral of -p n over the true plate, p sampled in the fluid and extrapolated to the wall	0.886	0.886	1.254
discrete divergence theorem over the marked staircase, p taken from the adjacent fluid cell	1.016	0.907	1.362
discrete divergence theorem over the marked staircase, p taken from the body cell (what the ghost fill holds)	-0.576	-0.270	0.658
MFC's own sum of -grad p dV over body cells (4th order, as in s_compute_ib_forces)	0.141	0.278	0.323

Read down that column. Using the fluid pressure just outside the body recovers the correct force to within about
9 percent, and it comes out normal to the plate as it must. Using the pressure that the ghost-cell reconstruction
has written into the body gives a force of the wrong sign. MFC's volume integral is a sum over exactly those
cells, so it inherits the error.

Two things this rules out:

• It is not the stencil width. Repeating the volume integral with a 3-point second-order stencil instead of
  the 5-point fourth-order one gives 0.393 instead of 0.323, against a true 1.254. Narrowing the stencil barely
  helps, so incomplete telescoping across a thin body is not the cause. (The body is 5 cells across along both x
  and z at this incidence, measured directly.)
• It is not the flow. The fluid-side pressures around the body give the right answer.

Why the interior pressure comes out inverted

I can measure the symptom but have not read enough of the ghost-point machinery to be certain of the cause. The
shape of it suggests the image point: a ghost cell inside a thin plate is reflected across its nearest surface by
twice its wall distance, and for a body only a few cells thick that image can land on, or beyond, the opposite
face. The ghost then takes its value from the far side of the plate, which inverts the pressure jump exactly as
observed. That would also explain why the error shrinks but does not vanish as the section is refined: fewer
cells are close enough to the far face, but the ones near mid-thickness still are.

A fix that the measurement supports

Compute the force from the fluid side rather than the body side: for every marked cell, for each of its faces
whose neighbour is fluid, accumulate -p_fluid n dA (and the corresponding viscous traction). That is the same
discrete divergence theorem the current sum is trying to evaluate, but it reads the pressure where the pressure
is trustworthy. On this case it gives 1.362 against a true 1.254, where the present code gives 0.323.

It is also cheaper: one pass over marked cells testing six neighbours, with no finite-difference stencil and no
fd_order dependence.

I am happy to prepare that as a PR if the approach looks right to you. The reason I have not simply pushed one is
that it changes the reported force on every existing immersed-boundary case, and the golden files that go with
them, so it seems worth agreeing on before doing.

Reproducing

tools/dissect_ib_force.py in the case repository does the decomposition above from a post-processed field; the
three-line summary is the table. examples/3D_ibm_pitchup_plate (PR #1847) is the case, and its C_D / C_L
ratio after the ramp is a reproducer needing no reference data.