2nd-order FDM discretization of this equation system on a regular mesh (with rectangular elements) yields a 9-noded molecule in variable 'u' with the usual W, E, S, N, and C (center) coeffs, but also NW, NE, SW, and SE.

My question is how to correctly apply a zero flux bc, especially when acute or obtuse corners are present on the mesh boundaries of an irregular-shaped 2D domain?

(1). I guess there is the missing part. That is the mesh. (2). By using the coordinate transformation, you will have the governing equation written in the general coordinates, which is yet unknown. (3). But if you then create a mesh which is normal to the boundary, then the zero flux boundary condition is just: t(Xi_at wall, Eta, Zeta)= t(Xi_at wall - 1, Eta, Zeta), if Xi is the coordinate normal to the local wall. (4). If you are serious about the zero flux condition, then the numerical grid generation is a good one to try. (5). By the way, I normally by-pass the cooking step, and buy the Chinese food (cooked) at the local mall's PANDA EXPRESS fast food place. (6). So, zero flux means the normal gradient is zero, that is you need to define the local wall point and a point on the line normal to the wall point. That means a local mesh line must be normal to the wall. When you have this local behavior of the mesh, the zero flux condition is very simple.

(1). If on the other hand, you insist on using the regular coordinates (Cartesian or the like), then you can find a point of the mesh near the wall first, then draw a normal to the wall, find the wall point, then assign the value on the wall equal to the near the wall point. (2). After you have these wall point values updated, then use interpolation to compute the values at other wall points if needed.